1. Loading data

1.1. phyloseq obejct

1.2. qPCR data (controls)

2. QC

QC1. How many samples failed sequencing

QC2. How were changes in read stats and host DNA proportion?

QC3. How were the extraction controls

QC4. Prevalence / abundance filtering - red flag

3. Analysis

A0. Calculation of alpha-diversity indices

A1. Host DNA, bacterial DNA and % host

A2. Modeling of sequencing results

A3. Taxa alpha diversity

A4. Mediation analysis on species richness

A5. Taxa beta diversity

A6. DA for taxa

A7. LM of function alpha diversity

A8. permanova of function alpha diversity

A9. DA for function

A10. Rarefaction analysis

Data inputs

Meta data

• qPCR - bacteria

• qPCR - human

• qPCR host %

• Raw reads

• final reads

• sequencing host %

• library prep failure status

• Raw reads

• subject_id

• treatment

• sample_type

• subject_id

Sequencing result

• samples

• controls

QC table by treated (binary)

Changes in matrices were observed

#sequencing result by sample type and control (1/0)


sample_data %>% data.frame() %>% 
        dplyr::group_by(sample_type, treated) %>% 
        dplyr::summarise(N = n(),
                  `Raw reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Raw_reads/10000000),2),nsmall = 2, big.mark = ","), " [", format(round(quantile(Raw_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Raw_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
                  `Host reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Host_mapped/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Host_mapped/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Host_mapped/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
                  `Host reads proportion<br>(median [IQR])<br>[%]` = paste(format(round(median(sequencing_host_prop * 100),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(sequencing_host_prop * 100, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(sequencing_host_prop * 100, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
                  `Final reads<br>(median [IQR])<br>[reads x 10<sup>7</sup>]` = paste(format(round(median(Final_reads/10000000),2), nsmall = 2, big.mark = ","), " [", format(round(quantile(Final_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Final_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), "]", sep = ""),
        ) %>%
        dplyr::rename(`Sample type` = sample_type, Treated = treated) %>%
        data.frame(check.names = F) %>% mutate_all(linebreak) %>% kbl(format = "html", escape = F) %>% kable_styling(full_width = 0, html_font = "serif")

QC table by treatment methods

Changes were sample type * treatment specific

table1 <- sample_data %>% data.frame() %>% 
        dplyr::filter(sample_type %in% c("Sputum", "Nasal", "BAL")) %>% 
        group_by (sample_type, treatment) %>%
        dplyr::summarise(N = n(),
              `Total DNA <br>ng/µL` = paste(format(round(median(DNA_host_ng_uL + DNA_bac_ng_uL),2), nsmall = 2, big.mark = ","), "<br>(", format(round(quantile(DNA_host_ng_uL + DNA_bac_ng_uL, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(DNA_host_ng_uL + DNA_bac_ng_uL, 0.75),2), nsmall = 2, big.mark = ","), ")", sep = ""),
              `Library or sequencing failure<br>N (%)` = paste(sum(lib_failed), "<br>(", sum(lib_failed) / n() * 100, " %)", sep = ""),
              `Clean reads<br>reads x 10<sup>7</sup>` = paste(format(round(median(Reads_after_trim/10000000),2), nsmall = 2, big.mark = ","), "<br>(", format(round(quantile(Reads_after_trim/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Reads_after_trim/10000000, 0.75), 1), nsmall = 1, big.mark = ","), ")", sep = ""),
              `Host reads %` = paste(format(round(median(sequencing_host_prop) * 100, 1), nsmall = 1, big.mark = ","), "<br>(", format(round(quantile(sequencing_host_prop, 0.25) * 100, 1), nsmall = 1, big.mark = ","), ", ", format(round(quantile(sequencing_host_prop, 0.75) * 100, 1), nsmall = 1, big.mark = ","), ")", sep = ""),
              `Final reads<br>reads x 10<sup>7</sup>` = paste(format(round(median(Final_reads/10000000),2), nsmall = 2, big.mark = ","), "<br>(", format(round(quantile(Final_reads/10000000, 0.25),2), nsmall = 2, big.mark = ","), ", ", format(round(quantile(Final_reads/10000000, 0.75),2), nsmall = 2, big.mark = ","), ")", sep = "")
        ) %>% data.frame(check.names = F) %>% 
        arrange(sample_type, treatment) %>%
        dplyr::rename(`Sample type` = sample_type, Treatment = treatment) %>%
        mutate_all(linebreak) %>% kbl(format = "html", escape = F) %>% kable_styling(full_width = 0, html_font = "serif")

table1 

Figure of reads by treatment (z-score)

Changes were sample type * treatment specific

# Summary figures - facet and z-score -------------------------------------

sample_data %>% 
        subset(., !is.na(.$subject_id)) %>%
        data.frame() %>%
        gather(feature, value, Raw_reads:sequencing_host_prop) %>%
        group_by(feature, sample_type) %>% 
        subset(., .$feature %in% c("Raw_reads", "Host_mapped", "Final_reads", "sequencing_host_prop")) %>%
        mutate(z_score = scale(value),
               feature = factor(feature, levels = c("Raw_reads", "Host_mapped", "Final_reads", "sequencing_host_prop"), labels = c("Raw reads", "Host mapped", "Final reads", "Host %"))) %>%
        ggplot(aes(x = treatment, y = z_score, fill = treatment)) +
                geom_boxplot(lwd = 0.2) +
                guides(fill=guide_legend(title="Treatment", nrow = 1)) +
                facet_grid(sample_type~feature) +
                xlab("Treatment") +
                ylab("Z score") +
                theme_classic(base_family = "serif", base_size = 14) +
                guides( x =  guide_axis(angle = 90)) + 
                theme(legend.position = "top") +
                scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment") #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6

Results:

2.1. There were no differences in raw reads.

2.2. However, final reads increased after some treatment, and host DNA proportion decreased

Species richness of controls

Some possible contaminants were identified in extraction controls

#Loading theoretical mock community
zymo_mock <- read_excel("/Users/minsikkim/Dropbox (Partners HealthCare)/@minsik/project_sicas2/data_raw/DAR_20210929_zymo_mock_data.xlsx") %>%
        data.frame(row.names = T) %>% rename(mock_theoretical = Mock) %>% mutate(mock_theoretical = mock_theoretical/100) %>%
        merge_phyloseq(otu_table(., taxa_are_rows = T), tax_table(phyloseq$phyloseq_count))

phyloseq_mock <- rbind(c("mock_theoretical", "Mock theoretical", "-")) %>% data.frame() %>%
        column_to_rownames(var = "X1") %>% rename(sample_type = X2, treatment = X3) %>% #making sample_data of mock community
        merge_phyloseq(sample_data(.), zymo_mock)


phyloseq_control_rel <- subset_samples(phyloseq$phyloseq_rel, sample_type == "Mock" | sample_type == "Neg.") #adding data of controls
sample_data(phyloseq_control_rel)$treatment <- sample_data(phyloseq_control_rel)$treatment %>% as.character()
sample_data(phyloseq_control_rel)$sample_type <- sample_data(phyloseq_control_rel)$sample_type %>% as.character()
phyloseq_control_rel <- merge_phyloseq(phyloseq_control_rel, phyloseq_mock) 


#Species richness of each control groups
sample_data(phyloseq_control_rel)$S.obs <- rowSums(t(otu_table(phyloseq_control_rel)) != 0)
sample_data(phyloseq_control_rel)$sample_type <- 
        factor(sample_data(phyloseq_control_rel)$sample_type, levels = c("Mock theoretical", "Mock", "Neg."))
sample_data(phyloseq_control_rel)$teratment <-
        factor(sample_data(phyloseq_control_rel)$treatment, levels = c("-", "Control", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp"))
phyloseq_control_rel %>%
        sample_data() %>%
        #mutate(sample_type = factor(sample_type, levels = c("Mock", "Neg.")),
        #       treatment = factor(treatment, levels = c("Theoretical", "Control", "Benzonase", "Host zero", "Molysis", "QIAamp"))) %>%
        group_by(sample_type, treatment) %>%
        summarise(Mean = mean(S.obs),
                  SD = sd(S.obs),) %>%
        kbl(format = "html", caption = "Species richness of controls") %>% 
        kable_styling(full_width = 0, html_font = "serif")

Bar plot of controls

Some possible contaminants were identified in extraction controls

Some changes visible at postive control….

#Manipulating phyloseq - only top 10 
tax_table(phyloseq_control_rel) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq_control_rel) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 8] <- .[, 8] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq_control_rel
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        plot_bar(., fill="species20") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Bar plot of control data") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ sample_type, scales= "free_x", nrow=1)

#there could be opportunistic pathogens...

Bar plot of controls (Positive)

Some possible contaminants were identified in extraction controls

Gram negatives were fragile to depletion method at postivive control

#Manipulating phyloseq - only top 10 

phyloseq_control_rel %>% 
        subset_samples(., sample_type == "Mock") %>% 
        tax_table() %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(subset_samples(phyloseq_control_rel,sample_type == "Mock" & S.obs != 0)) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 8] <- .[, 8] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- subset_samples(phyloseq_control_rel,sample_type == "Mock" & S.obs != 0)
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        plot_bar(., fill="species20") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Postive controls") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Control", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

#there could be opportunistic pathogens...

Issue 1 - gram negative

Benzonse and Host zero depleted all the gram negative strains

Others decrased gram negatives a lot, but not became zero

#Making function without borders
my_plot_bar = function (physeq, x = "Sample", y = "Abundance", fill = NULL, title = NULL, 
                        facet_grid = NULL) {
    mdf = psmelt(physeq)
    p = ggplot(mdf, aes_string(x = x, y = y, fill = fill))
    p = p + geom_bar(stat = "identity")
    p = p + theme(axis.text.x = element_text(angle = -90, hjust = 0))
    if (!is.null(facet_grid)) {
        p <- p + facet_grid(facet_grid)
    }
    if (!is.null(title)) {
        p <- p + ggtitle(title)
    }
    return(p)
}
# gram stain data
phyloseq_control_rel %>% 
        subset_samples(., sample_type == "Mock") %>% 
        tax_table() %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(subset_samples(phyloseq_control_rel,sample_type == "Mock" & S.obs != 0)) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 8] <- .[, 8] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- subset_samples(phyloseq_control_rel,sample_type == "Mock" & S.obs != 0)
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        my_plot_bar(., fill="gram_stain") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Gram stain in Zymo mock") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Gram-stain")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

#there could be opportunistic pathogens...
                                
#Manipulating phyloseq - only top 10 

sample_data(phyloseq_control_rel) <- cbind(phyloseq_control_rel %>%
              sample_data %>%
              data.frame(),
      phyloseq_control_rel %>% 
        otu_table %>%
        data.frame %>%
        subset(., rownames(.) %in% head(taxa_sums(subset_samples(phyloseq_control_rel,sample_type == "Mock" & S.obs != 0)) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)) %>%
        t()
)


sample_data(phyloseq_control_rel) %>%
        data.frame() %>%
        subset(., !is.na(.$Escherichia_coli)) %>% 
        group_by(sample_type, treatment) %>% 
        summarise(N = n(),
                  `<i>Escherichia coli</i><br>(median [IQR])` = paste(format(round(median(Escherichia_coli),4),nsmall = 4, big.mark = ","), " [", format(round(quantile(Escherichia_coli, 0.25),4), nsmall = 4, big.mark = ","), ", ", format(round(quantile(Escherichia_coli, 0.75),4), nsmall = 4, big.mark = ","), "]", sep = ""),
                  `<i>Pseudomonas aeruginosa</i><br>(median [IQR])` = paste(format(round(median(Pseudomonas_aeruginosa_group),4),nsmall = 4, big.mark = ","), " [", format(round(quantile(Pseudomonas_aeruginosa_group, 0.25),4), nsmall = 4, big.mark = ","), ", ", format(round(quantile(Pseudomonas_aeruginosa_group, 0.75),4), nsmall = 4, big.mark = ","), "]", sep = ""),
                  `<i>Salmonella enterica</i><br>(median [IQR])` = paste(format(round(median(Salmonella_enterica),4),nsmall = 4, big.mark = ","), " [", format(round(quantile(Salmonella_enterica, 0.25),4), nsmall = 4, big.mark = ","), ", ", format(round(quantile(Salmonella_enterica, 0.75),4), nsmall = 4, big.mark = ","), "]", sep = "")
        ) %>%
        rename(`Sample type` = sample_type) %>%
        data.frame(check.names = F) %>% mutate_all(linebreak) %>% kbl(format = "html", escape = F) %>% kable_styling(full_width = 0, html_font = "serif")

Issue 2 - Positive controls contaminants

Some possible contaminants were identified in most of samples

This could be 1) background contamination or 2) cross-contamination from kingfisher. Most of these are gram positives. Negative controls should be double-checked

#Manipulating phyloseq - only top 10 
phyloseq_control_rel_contam <- subset_taxa(phyloseq_control_rel , !(taxa_names(phyloseq_control_rel) %in% head(taxa_sums(subset_samples(phyloseq_control_rel,
                                                                                         sample_type == "Mock" & S.obs != 0)) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10))
)

phyloseq_control_rel_contam <- subset_taxa(phyloseq_control_rel_contam, taxa_sums(phyloseq_control_rel_contam) != 0)
phyloseq_control_rel_contam <- subset_samples(phyloseq_control_rel_contam, sample_type != "Neg." & S.obs != 0)


tax_table(phyloseq_control_rel_contam) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq_control_rel_contam) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 9] <- .[, 9] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq_control_rel_contam
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        my_plot_bar(., fill="species20") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Contaminants in Zymo mock") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Control", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

Negative controls

Contaminants of positive and negative control do not match

It seems without host DNA, gram-negatives are vulnerable to depletion methods.

These negative contaminants highly likely introduced after-depletion

#Manipulating phyloseq - only top 10 


tax_table(phyloseq_control_rel) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq_control_rel) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 9] <- .[, 9] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq_control_rel
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        subset_samples(., sample_type == "Neg.") %>%
        my_plot_bar(., fill="species20") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Barplot of neg. data") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

phyloseq_control_rel %>%
        subset_samples(., sample_type == "Neg.") %>%
        my_plot_bar(., fill="gram_stain") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Gram-stain of negative data") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

Samples - gram-stain?

Nasal swab had low gram positive/pattern consisted after depletion

BAL showed Similar gram - / + ratio

Sputum showed high decrease in gram negative bacteria

Freeze/thaw cycle could be associated

Currently no further analysis is possible

#Manipulating phyloseq - only top 10 

tax_table(phyloseq$phyloseq_rel) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq$phyloseq_rel) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 9] <- .[, 9] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq$phyloseq_rel
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        subset_samples(., sample_type == "Nasal") %>%
        my_plot_bar(., fill="gram_stain") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Gram stain of nasal samples") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

tax_table(phyloseq$phyloseq_rel) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq$phyloseq_rel) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 9] <- .[, 9] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq$phyloseq_rel
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        subset_samples(., sample_type == "BAL") %>%
        my_plot_bar(., fill="gram_stain") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Gram stain of BAL samples") +
        theme(legend.text = element_markdown()) +
        guides(fill=guide_legend(title="Top 10 species")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

tax_table(phyloseq$phyloseq_rel) %>%
        cbind(species20 = "[Others]") %>%
        {top20species <- head(taxa_sums(phyloseq$phyloseq_rel) %>%
                                data.frame %>%
                                arrange(-.) %>%
                                row.names(), 10)
   .[top20species, "species20"] <- as.character(.[top20species, "Species"])
   .[, 9] <- .[, 9] %>% gsub("s__", "", .) %>% gsub("_", " ", .) %>% paste("<i>", ., "</i>", sep = "")
   phyloseq_temp <- phyloseq$phyloseq_rel
   tax_table(phyloseq_temp) <- tax_table(.) 
   phyloseq_temp
  } %>%
        subset_samples(., sample_type == "Sputum") %>%
        my_plot_bar(., fill="gram_stain") + 
        ylab("Relative abundancne") +
        theme_classic(base_size = 11, base_family = "serif") +
        ggtitle("Gram stain of sputum samples") +
        theme(legend.text = element_markdown()) +
        scale_fill_brewer(type = "qual", palette = 6) +
        guides(fill=guide_legend(title="Gram-stain")) +
        facet_wrap (~ factor(treatment, levels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAamp")),
                    scales= "free_x", nrow=1)

Histogram of prelanence taxa

No prevalence or abundance filtering (each experimental group is 5% of total sample)

Before anlayzing, alpha diversity indices were calculated for all phyloseq objects

alpha_diversity <- function(data) {
        otu_table <- otu_table(data) %>% .[colSums(.) !=0]
        S.obs <- rowSums(t(otu_table) != 0)
        sample_data <- sample_data(data)
        data_evenness <- vegan::diversity(t(otu_table)) / log(vegan::specnumber(t(otu_table))) # calculate evenness index using vegan package
        data_shannon <- vegan::diversity(t(otu_table), index = "shannon") # calculate Shannon index using vegan package
        data_hill <- exp(data_shannon)                           # calculate Hills index
        
        data_dominance <- microbiome::dominance(otu_table, index = "all", rank = 1, aggregate = TRUE) # dominance (Berger-Parker index), etc.
        data_invsimpson <- vegan::diversity(t(otu_table), index = "invsimpson")                          # calculate Shannon index using vegan package
        alpha_diversity <- cbind(S.obs, data_shannon, data_hill, data_invsimpson, data_evenness,data_dominance) # combine all indices in one data table
        sample_data <- merge(data.frame(sample_data), alpha_diversity, by = 0, all = T) %>% column_to_rownames(var = "Row.names")
}

phyloseq_unfiltered <- phyloseq

sample_data(phyloseq_unfiltered$phyloseq_rel) <- sample_data(alpha_diversity(phyloseq_unfiltered$phyloseq_count))
sample_data(phyloseq_unfiltered$phyloseq_count) <- sample_data(alpha_diversity(phyloseq_unfiltered$phyloseq_count)) 
sample_data(phyloseq_unfiltered$phyloseq_path_rpk) <- sample_data(alpha_diversity(phyloseq_unfiltered$phyloseq_path_rpk))

#Calculation of sample prevalence, standard deviation, median abundance across all samples for all bugs and making into a table.
#
#•  In initial analysis we will not perform prevalence or abundance filtering (though we may consider this for secondary differential abundance analyses to manage p (features) > n (sample size) problem and issues with multiple hypothesis correction)


taxa_qc <- data.frame("species" =  otu_table(subset_samples(phyloseq_unfiltered$phyloseq_count,
                                                            S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum"))) %>%
                              t() %>% colnames(),
        "prevalence" = ifelse(subset_samples(phyloseq_unfiltered$phyloseq_count, S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum")) %>% otu_table() > 0, 1, 0) %>% t() %>% colSums(), #Prevalence of taxa
        "mean_rel_abd" = subset_samples(phyloseq_unfiltered$phyloseq_count, S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum")) %>%
                otu_table() %>%
                t() %>%
                colMeans(na.rm = T) #mean relativ abundacne 
)

function_qc <- data.frame("function" =  otu_table(subset_samples(phyloseq_unfiltered$phyloseq_path_rpk, S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum"))) %>% t() %>% colnames(),
        "prevalence" = ifelse(subset_samples(phyloseq_unfiltered$phyloseq_path_rpk, S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum")) %>% otu_table() > 0, 1, 0) %>% t() %>% colSums(), #Prevalence of taxa
        "mean_rpk" = subset_samples(phyloseq_unfiltered$phyloseq_path_rpk, S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum")) %>% otu_table() %>% t() %>% colMeans(na.rm = T) #mean relativ abundacne 
)


hist(log10(taxa_qc$prevalence), xlab = "log10(Taxa prevalence)", main = "Histogram of prevalence of taxa")

Histogram of mean abundance

hist(log10(taxa_qc$mean_rel_abd), xlab = "log10(Mean relative abundance)", main = "Histogram of mean relative abundance")

decontam - stratified by sample type

Decontam analysis result stratified by sample type

Decontam analysis result stratified by sample type

#Stratified by sample type
sample_data(phyloseq_unfiltered$phyloseq_rel)$is.neg <- grepl("Neg", sample_data(phyloseq_unfiltered$phyloseq_rel)$sample_type)

contaminant <- data.frame() 

#cat("decontam prevalence - all")
contaminant1 <- 
data.frame("prevalence", "all", fix.empty.names = F, 
phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0 & sample_type != "Mock") %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.1) %>% subset(.,.$contaminant) %>% row.names()
)

#cat("decontam prevalence - BAL")
contaminant2 <- 
data.frame("prevalence", "BAL", fix.empty.names = F, 
subset_samples(phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0), sample_type %in% c("BAL", "Neg.")) %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.1) %>% subset(.,.$contaminant) %>% row.names())

#cat("decontam prevalence - Nasal")
contaminant3 <- 
        data.frame("prevalence", "Nasal swab", fix.empty.names = F, 
subset_samples(phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0), sample_type %in% c("Nasal", "Neg.")) %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.1) %>% subset(.,.$contaminant) %>% row.names())

#cat("prevalence", "decontam prevalence - Sputum")
contaminant4 <- 
data.frame("prevalence", "Sputum", fix.empty.names = F,
subset_samples(phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0), sample_type %in% c("Sputum", "Neg.")) %>%
        isContaminant(., method="prevalence", neg = "is.neg", threshold = 0.1) %>% subset(.,.$contaminant) %>% row.names())

#cat("decontam frequency - All")
contaminant5 <- 
data.frame("frequency", "all", fix.empty.names = F,
phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0) %>%
        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant) %>% row.names())

#cat("decontam frequency - BAL")
#contaminant6 <- 
#data.frame("frequency", "BAL", fix.empty.names = F,
#subset_samples(phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0), sample_type %in% c("BAL")) %>%
#        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant) %>% row.names()
#)

#cat("decontam frequency - Nasal")
#contaminant7 <- 
#data.frame("frequency", "Nasal swab", fix.empty.names = F,
#subset_samples(phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0), sample_type %in% c("Nasal")) %>%
#        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant) %>% row.names()
#)

#cat("decontam frequency - Sputum")
contaminant8 <- 
data.frame("frequency", "Sputum", fix.empty.names = F,
subset_samples(phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0), sample_type %in% c("Sputum")) %>%
        isContaminant(method="frequency", conc="DNA_bac_well") %>% subset(.,.$contaminant) %>% row.names()
)
#cat("decontam combined - All")
contaminant9 <- 
data.frame("combined", "all", fix.empty.names = F,

phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0) %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant) %>% row.names())

#cat("decontam combined - BAL")
contaminant10 <- 
data.frame("combined", "BAL", fix.empty.names = F,
subset_samples(phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0), sample_type %in% c("BAL", "Neg.")) %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant) %>% row.names())

#cat("decontam combined - Nasal")
#contaminant11 <- 
#data.frame("combined", "Nasal swab", fix.empty.names = F,
#subset_samples(phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0), sample_type %in% c("Nasal", "Neg.")) %>%
#        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant) %>% row.names() 
#)

#cat("decontam combined - Sputum")
contaminant12 <- 
data.frame("combined", "Sputum", fix.empty.names = F,
subset_samples(phyloseq_unfiltered$phyloseq_rel %>% subset_samples(S.obs != 0), sample_type %in% c("Sputum", "Neg.")) %>%
        isContaminant(method="combined", neg="is.neg", conc = "DNA_bac_well") %>% subset(.,.$contaminant) %>% row.names())

#Some contaminant table was not included as they were empty. (Contaminant 6, 7, 11)
contaminant <- rbind(contaminant1, contaminant2, contaminant3, contaminant4, contaminant5, contaminant8, contaminant9, contaminant10, contaminant12)
names(contaminant) <- c("method", "sample_type", "contaminant")
cat("List of contaminant (prevalence method) - from stratified analysis\n\n")

## List of contaminant (prevalence method) - from stratified analysis

contaminant %>% subset(., .$method == "prevalence" & .$sample_type != "all") %>% .$contaminant %>% unique()

##  [1] "Cutibacterium_acnes"           "Thermoleophilum_album"        
##  [3] "Paludisphaera_borealis"        "Escherichia_coli"             
##  [5] "Bacillus_intestinalis"         "Listeria_monocytogenes"       
##  [7] "Microbacterium_ginsengisoli"   "Microbacterium_laevaniformans"
##  [9] "Arthrobacter_agilis"           "Brevundimonas_diminuta"       
## [11] "Bosea_vaviloviae"              "Acetobacter_aceti"            
## [13] "Alcaligenes_faecalis"          "Pseudomonas_formosensis"      
## [15] "Enterococcus_faecalis"         "Lactobacillus_fermentum"

cat("Prevalence decontam table\n\n")

## Prevalence decontam table

contaminant %>% subset(., .$method == "prevalence") %>% .$sample_type %>% table

## .
##        all        BAL Nasal swab     Sputum 
##         16         14         12         12

cat("Frequency decontam table\n\n")

## Frequency decontam table

contaminant %>% subset(., .$method == "frequency") %>% .$sample_type %>% table

## .
##    all Sputum 
##      3      3

cat("Combined decontam table\n\n")

## Combined decontam table

contaminant %>% subset(., .$method == "combined") %>% .$sample_type %>% table

## .
##    all    BAL Sputum 
##      6      1      2

matrix(nrow=3,ncol=4) %>% data.frame() %>% rename(" " = X1, "BAL" = X2, "Nasal swab" = X3, "Sputum" = X4) %>%
        rownames_to_column("x") %>% mutate(x = c("Prevalence method", "Frequency method", "Combined"),
                                           " " = c(16, 3, 6),
                                           "BAL" = c(14, 0, 1),
                                           "Nasal swab" = c(12, 0, 0),
                                           "Sputum" = c(12, 3, 2)
                                           ) %>% column_to_rownames("x") %>%
        kbl(format = "html") %>%
        add_header_above(c(" " = 1, "Non-stratified" = 1, "Stratified" = 3)) %>% 
        kable_styling(full_width = 0, html_font = "serif") 

Red flags

Taxa with low prevalences were red-flagged

Prevalence/abundance filtering was employed

taxa_qc <- data.frame("species" =  otu_table(subset_samples(phyloseq_unfiltered$phyloseq_count,
                                                            S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum"))) %>%
                              t() %>% colnames(),
        "prevalence" = ifelse(subset_samples(phyloseq_unfiltered$phyloseq_count, S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum")) %>% otu_table() > 0, 1, 0) %>% t() %>% colSums(), #Prevalence of taxa
        "mean_rel_abd" = subset_samples(phyloseq_unfiltered$phyloseq_count, S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum")) %>%
                otu_table() %>%
                t() %>%
                colMeans(na.rm = T) #mean relativ abundacne 
)

function_qc <- data.frame("function" =  otu_table(subset_samples(phyloseq_unfiltered$phyloseq_path_rpk, S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum"))) %>% t() %>% colnames(),
        "prevalence" = ifelse(subset_samples(phyloseq_unfiltered$phyloseq_path_rpk, S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum")) %>% otu_table() > 0, 1, 0) %>% t() %>% colSums(), #Prevalence of taxa
        "mean_rpk" = subset_samples(phyloseq_unfiltered$phyloseq_path_rpk, S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum")) %>% otu_table() %>% t() %>% colMeans(na.rm = T) #mean relativ abundacne 
)

red_flag_taxa <- data.frame(species = taxa_qc$species,
                            red_flag_prev_abd = ifelse(taxa_qc$prevalence < otu_table(phyloseq_unfiltered$phyloseq_rel) %>%
                                                               t %>% rownames() %>%
                                length * 0.05 & taxa_qc$mean_rel_abd < quantile(taxa_qc$mean_rel_abd, 0.75), 1,0)) %>%
        mutate(red_flag_decontam_prev = species %in% (contaminant %>%
                                    subset(., .$method == "prevalence" & .$sample_type != "all") %>%
                                    .$contaminant %>% unique()))


red_flag_function <- data.frame(function. = function_qc$function., red_flag_prev_abd = ifelse(function_qc$prevalence < otu_table(phyloseq$phyloseq_path_rpk) %>% t %>% rownames() %>% length * 0.05 & function_qc$mean_rpk < quantile(function_qc$mean_rpk, 0.75), 1, 0))

#phyloseq of table 1
phyloseq$phyloseq_count_prev <- prune_taxa(subset(red_flag_taxa,
                                           red_flag_taxa$red_flag_decontam_prev != 1)$species,
                                    phyloseq$phyloseq_count)
#phyloseq for analysis
phyloseq$phyloseq_rel <- prune_taxa(subset(red_flag_taxa,
                                           red_flag_taxa$red_flag_prev_abd != 1 & !red_flag_taxa$red_flag_decontam_prev)$species,
                                    phyloseq$phyloseq_rel)
phyloseq$phyloseq_count <- prune_taxa(subset(red_flag_taxa,
                                           red_flag_taxa$red_flag_prev_abd != 1 & !red_flag_taxa$red_flag_decontam_prev)$species,
                                    phyloseq$phyloseq_count)
phyloseq$phyloseq_path_rpk <- prune_taxa(subset(red_flag_taxa, red_flag_function$red_flag_prev_abd != 1)$function., phyloseq$phyloseq_path_rpk)

#phyloseq$tree_phyloseq_count <- prune_taxa(subset(red_flag_taxa,
                                           #red_flag_taxa$red_flag_prev_abd != 1 & !red_flag_taxa$red_flag_decontam_prev)$species,
                                    #phyloseq$tree_phyloseq_count)

#phyloseq$tree_phyloseq_rel <- prune_taxa(subset(red_flag_taxa,
                                           #red_flag_taxa$red_flag_prev_abd != 1 & !red_flag_taxa$red_flag_decontam_prev)$species,
                                    #phyloseq$tree_phyloseq_rel)

Before anlayzing, alpha diversity indices were calculated for all filtred phyloseq objects

#Calculation of alpha diversity indices for filtered samples

alpha_diversity <- function(data) {
        otu_table <- otu_table(data) %>% .[colSums(.) !=0]
        S.obs <- rowSums(t(otu_table) != 0)
        sample_data <- sample_data(data)
        data_evenness <- vegan::diversity(t(otu_table)) / log(vegan::specnumber(t(otu_table))) # calculate evenness index using vegan package
        data_shannon <- vegan::diversity(t(otu_table), index = "shannon") # calculate Shannon index using vegan package
        data_hill <- exp(data_shannon)                           # calculate Hills index
        data_dominance <- microbiome::dominance(otu_table, index = "all", rank = 1, aggregate = TRUE) # dominance (Berger-Parker index), etc.
        data_invsimpson <- vegan::diversity(t(otu_table), index = "invsimpson")                          # calculate Shannon index using vegan package
        alpha_diversity <- cbind(S.obs, data_shannon, data_hill, data_invsimpson, data_evenness,data_dominance) # combine all indices in one data table
        sample_data <- merge(data.frame(sample_data), alpha_diversity, by = 0, all = T) %>% column_to_rownames(var = "Row.names")
}
sample_data(phyloseq$phyloseq_count_prev) <- sample_data(alpha_diversity(phyloseq$phyloseq_count_prev))
sample_data(phyloseq$phyloseq_rel) <- sample_data(alpha_diversity(phyloseq$phyloseq_count))
sample_data(phyloseq$phyloseq_count) <- sample_data(alpha_diversity(phyloseq$phyloseq_count)) 
sample_data(phyloseq$phyloseq_path_rpk) <- sample_data(alpha_diversity(phyloseq$phyloseq_path_rpk))  
#sample_data(phyloseq$tree_phyloseq_count) <- sample_data(alpha_diversity(phyloseq$tree_phyloseq_count))
#phyloseq$tree_phyloseq_rel <- transform_sample_counts(phyloseq$tree_phyloseq_count, function (x) {x/sum(x)})
sample_data <- sample_data(phyloseq$phyloseq_count)

qPCR result

#2A: Change in total DNA (qPCR)

f2a <- ggplot(sample_data, aes(x = sample_type, y = log10(DNA_host_nondil + DNA_bac_nondil))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        ylab("log<sub>10</sub>(qPCR total DNA)<br>(ng/μL)") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif") + 
        labs(tag = "A") +
        #scale_x_discrete(label = c( "Mock", "Neg.", "BAL", "Nasal", "Sputum")) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))


#2B: Change in human DNA (qPCR)
f2b <- ggplot(sample_data, aes(x = sample_type, y = log10(DNA_host_nondil))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3")) +
        ylab("log<sub>10</sub>(qPCR host DNA)<br>(ng/μL)") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif")+ 
        labs(tag = "B") +
        #scale_x_discrete(label = c("BAL", "Nasal", "Sputum", "Mock", "Neg.")) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))
#2C: Change in 16S DNA (qPCR)
f2c <- ggplot(sample_data, aes(x = sample_type, y = log10(DNA_bac_nondil))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3")) +
        ylab("log<sub>10</sub>(qPCR bacterial DNA)<br>(ng/μL)") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif")+ 
        labs(tag = "C") +
        #scale_x_discrete(label = c("BAL", "Nasal", "Sputum", "Mock", "Neg.")) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))

#2D. Change in % host (qPCR)
f2d <- ggplot(sample_data, aes(x = sample_type, y = host_proportion)) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3")) +
        ylab("Host DNA ratio") +
        xlab("Sample type") +
        theme_classic (base_size = 12, base_family = "serif") + 
        labs(tag = "D") +
        #scale_x_discrete(label = c("BAL", "Nasal", "Sputum", "Mock", "Neg.")) +
        theme(plot.tag = element_text(size = 15), axis.title.y = element_markdown()) +              # Plot title size
        guides(fill = guide_legend(nrow = 1, title = "Treatment"))

#output for markdown
ggarrange(f2a, f2b, f2c, f2d, common.legend = T , align = "hv")

Figure 2. qPCR result of host depletion study. A. Total DNA B. Host DNA C. Bacterial DNA D. Host %

Sequencing result

f3a <- ggplot(sample_data, aes(x = sample_type, y = log10(Raw_reads))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        theme_classic (base_size = 12, base_family = "serif") + 
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("Control","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        scale_x_discrete(name ="Sample type") +
        theme(axis.title.y = element_markdown(),
              plot.tag = element_text(size = 15)) +
        ylab("log<sub>10</sub>(raw reads)") +
        labs(tag = "A") +
        guides(fill = guide_legend(nrow = 1))

#   - Host_mapped


f3b <- ggplot(sample_data, aes(x = sample_type, y = log10(Host_mapped))) +
        theme_classic (base_size = 12, base_family = "serif")+ 
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("Control","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        scale_x_discrete(name ="Sample type")+
        theme(axis.title.y = element_markdown(),
              plot.tag = element_text(size = 15)) +
        ylab("log<sub>10</sub>(host reads)") +
        labs(tag = "B") +
        guides(fill = guide_legend(nrow = 1))


#   - % Host (we have used Host_mapped/Raw_reads in prior papers)

#   - Final_reads

f3c <- ggplot(sample_data, aes(x = sample_type, y = log10(Final_reads))) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("Control","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        scale_x_discrete(name ="Sample type")+
        ylab("log<sub>10</sub>(final reads)") +
        theme_classic (base_size = 12, base_family = "serif") + 
        theme(axis.title.y = element_markdown(),
              plot.tag = element_markdown(size = 15)) +
        labs(tag = "C") +
        guides(fill = guide_legend(nrow = 1))


#   - % Host (we have used Host_mapped/Raw_reads in prior papers)
f3d <- ggplot(sample_data, aes(x = sample_type, y = sequencing_host_prop)) +
        geom_boxplot(aes(fill = treatment), lwd = 0.2) +
        theme_classic (base_size = 12, base_family = "serif")+ 
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("Control","lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        scale_x_discrete(name ="Sample type")+
        theme(axis.title.y = element_markdown(),
              plot.tag = element_text(size = 15)) +
        ylab("Host ratio by sequencing") +
        labs(tag = "D") +
        guides(fill = guide_legend(nrow = 1))



ggarrange(f3a, f3b, f3c, f3d, common.legend = T, align = "hv")

Figure 3. Sequencing result of host depletion study. A. Total DNA B. Host DNA C. Bacterial DNA D. Host %

Library failure - ANOVA

Some samples failed in library prep. What type of sample were fragile to treatments?

glm ( library fail ~ sample_type + treatment + sample_type * treatment + subject_id )

glmer(lib_failed ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        Anova %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>Chisq)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>% column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Library failure

glm ( sequencing fail ~ sample_type + treatment + sample_type * treatment + subject_id )

Nasals were fragile to lyPMA and Molysis

glmer(lib_failed ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`t value`) > 2 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("sample_type|treatment", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>% 
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Sequencing failure

Modeling of sequencing failure were not available due to low number of cases.

BAL079 - control & lyPMA failed sequencing.

sample_data(phyloseq$phyloseq_count) %>% data.frame %>% mutate(sequencing_fail = (S.obs == 0)) %>%
glmer(sequencing_fail ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = .) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`t value`) > 2 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("sample_type|treatment", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>% 
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

log10(Final reads) - ANOVA

Which methods was effective in increasing the final reads?

Interaction term was significant

lmer(log10(Final_reads) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        anova %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>% column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

log10(Final reads)

Which methods was effective in increasing the final reads?

lmer( log10(Final reads) vs sample_type + treatment + sample_type * treatment + subject_id )

Except lyPMA, every methods increased final reads

lmer(log10(Final_reads) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("sample_type|treatment", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>% 
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Final reads x 10^6 - ANOVA

Which methods was effective in increasing the final reads?

Interaction term was significant

lmer(Final_reads/1000000 ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        anova %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>% column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Final reads x 1000000

Which methods was effective in increasing the final reads?

lmer( log10(Final reads) vs sample_type + treatment + sample_type * treatment + subject_id )

Except lyPMA, every methods increased final reads

lmer(Final_reads/1000000 ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("sample_type|treatment", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>% 
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Host ratio ANOVA

Which methods was effective in lowering host %

Interaction term was significant

lmer(sequencing_host_prop ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        anova %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>% column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Host ratio

Which methods was effective in lowering host %

lmer( Host DNA ratio vs sample_type + treatment + sample_type * treatment + (1|subject_id) )

Host zero was effect to to all types. Molysis was effective to Nasal and sputum. QIAamp was effective for Nasal only.

lmer(sequencing_host_prop ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*", .default = " ")) %>%
        rownames_to_column(var = "x") %>% mutate(x = gsub("sample_type|treatment", "", x))  %>% mutate(x = gsub(":", " * ", x)) %>% 
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Gram negatives - ANOVA

Which methods was changed gram-strain ratio?

Square root transformation was required

Interaction term was significant

hist(sample_data %>% data.frame %>% .$gram_neg_prop)

hist(sample_data %>% data.frame %>% .$gram_neg_prop %>% sqrt())

lmer(sqrt(gram_neg_prop) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        anova %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*", .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>% column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Gram negatives

Which method biased gram positive-negative ratio

lmer( Gram-negative ratio vs sample_type + treatment + sample_type * treatment + (1|subject_id) )

Some treatment (commercial) changed gram negative proportion

lmer(sqrt(gram_neg_prop) ~ sample_type + treatment + sample_type * treatment + (1|subject_id), data = sample_data %>% data.frame) %>%
        summary %>% .$coefficients %>% data.frame(check.names = F) %>%
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*", .default = " ")) %>%
        rownames_to_column(var = "x") %>% mutate(x = gsub("sample_type|treatment", "", x))  %>% mutate(x = gsub(":", " * ", x)) %>% 
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>% kable_styling(full_width = 0, html_font = "serif")

Species richness (all samples & interaction term) - ANOVA

Interaction term was significant

sample_data <- sample_data(phyloseq$phyloseq_count) %>% data.frame(check.names = F) %>% subset(., !is.nan(.$simpson))

lmer_sob <- lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = sample_data)
lmer_sob %>% 
        anova() %>% 
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Species richness (all samples & interaction term)

Increase at sputum was at every treatment Postive and negative control showed no changes

lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = sample_data) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Species richness - stratified (Pos + Neg)

No treatment increased species richenss - after adjusting sequencing depth. With mock community except lyPMA, treatments showed they even reduced the possible contaminants. Need to observe alpha diversity of positive controls

lm(S.obs ~ sample_type * treatment + log10 (Final_reads), data = subset(sample_data, sample_data$sample_type == "Neg." | sample_data$sample_type == "Mock" )) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Species richness - stratified (NS + Pos + Neg)

Molysis and host zero may incrased speciess richness of Nasal Data include nasal swab, positive depletion, and negative depletion

lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Nasal" | sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Species richness (BAL + Pos + Neg)

No changes observed

lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "BAL"  | sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Species richness (sputum + Pos + Neg)

Benzonase may incrased speciess richness of sputum

lmer(S.obs ~ sample_type * treatment + log10 (Final_reads) + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "Sputum"  | sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inv Simp - ANOVA

Final reads did not affect inverse Simpson

lmer_invsimpson <- lmer(data_invsimpson ~ sample_type * treatment + log10(Final_reads) + (1|subject_id), data = sample_data)

lmer_invsimpson %>% 
        anova() %>% 
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

lmer_invsimpson <- lmer(data_invsimpson ~ sample_type * treatment + (1|subject_id), data = sample_data)

Simpson (all samples & interaction term)

Sputum after treatment showed differences - stratified analysis is required

#Simpson

lmer_invsimpson %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inverse Simpson - stratified (Controls)

Inverse Simpson ~ sample_type + log10 (Final_reads) + (1|original_sample)

Mock community treated with lyPMA showed cahnges in alpha diveresity

lmer(data_invsimpson ~ sample_type * treatment + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inverse Simpson - stratified (NS + Pos + Neg)

Inverse Simpson ~ sample_type + log10 (Final_reads) + (1|original_sample)

Mock community treated with lyPMA only showed changes in alpha diveresity.

lmer(data_invsimpson ~ sample_type * treatment + (1|subject_id), data = subset(sample_data, sample_data$sample_type == "Nasal" | sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inverse Simpson - stratified (BAL + Pos + Neg)

Nothing changed in BAL

lmer(data_invsimpson ~ sample_type * treatment + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "BAL" |sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

Inverse Simpson - stratified (spt + Pos + Neg)

Sputum changed after some treatment - but their changes were not treatment global.

lmer(data_invsimpson ~ sample_type * treatment + log10(Final_reads) + + (1|original_sample), data = subset(sample_data, sample_data$sample_type == "Sputum" |sample_data$sample_type == "Mock" | sample_data$sample_type == "Neg.")) %>% 
        summary() %>%
        .$coefficients %>%
        data.frame(check.names = F) %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rownames_to_column(var = "x") %>% mutate(x = gsub("treatment|sample_type", "", x)) %>% mutate(x = gsub(":", " * ", x)) %>%
        column_to_rownames(var = "x") %>% 
        kbl(format = "html") %>%
        kable_styling(full_width = 0, html_font = "serif")

PCoA based on Horn-Morichita (all samples)

Jaccard dissimilarities (presenece and absence) showed BAL and Mock communities are distant. Some samples may have some overlaps

Changed to Horn-Morisita. It is less sensitive to changes due to sequencing depth

phyloseq_rel_nz <- transform_sample_counts(phyloseq$phyloseq_path_rpk, function(x) {x/sum(x)}) %>%
        subset_samples(S.obs != 0 & sample_type %in% c("Mock", "BAL", "Nasal", "Sputum"))

horn_perm_all <- vegan::adonis2(distance(phyloseq_rel_nz, method="horn") ~ sample_type * treatment + subject_id + log10(Final_reads),
                                  data = phyloseq_rel_nz %>% sample_data %>% data.frame(check.names = F),
                                  permutations = 10000)

horn_perm_ns <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "Nasal"), method="horn") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads),
                               data = subset_samples(phyloseq_rel_nz, sample_type == "Nasal") %>%
                                       sample_data %>% data.frame(check.names = F),
                               strata = subset_samples(phyloseq_rel_nz, sample_type == "Nasal") %>% 
                                       sample_data %>% data.frame(check.names = F) %>% .$subject_id, permutations = 10000)

horn_perm_bal  <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "BAL"), method="horn") ~  lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads),
                                 data = subset_samples(phyloseq_rel_nz, sample_type == "BAL") %>% sample_data %>% data.frame(check.names = F),
                                 strata = subset_samples(phyloseq_rel_nz, sample_type == "BAL") %>%
                                         sample_data %>% data.frame(check.names = F) %>% .$subject_id,
                                  permutations = 10000)

horn_perm_spt <- vegan::adonis2(distance(subset_samples(phyloseq_rel_nz, sample_type == "Sputum"), method="horn") ~ lypma + benzonase + host_zero + molysis + qiaamp + log10(Final_reads),
                                data = subset_samples(phyloseq_rel_nz, sample_type == "Sputum") %>% sample_data %>% data.frame(check.names = F),
                                strata = subset_samples(phyloseq_rel_nz, sample_type == "Sputum")
                                %>% sample_data %>% data.frame(check.names = F) %>% .$subject_id,
                                  permutations = 10000)



horn_perm_all %>% data.frame(check.names = F) %>% rownames_to_column("row.names") %>% 
        mutate(row.names = case_when(row.names == "sample_type" ~ 'Sample type',
                                     row.names == "treatment" ~ 'Treatment',
                                     row.names == "subject_id" ~ 'Subject',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "sample_type:treatment" ~ 'Sample type * Treatment',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
        mutate(` ` = case_when(abs(`Pr(>F)`) < 0.001 ~ "***",
                                            abs(`Pr(>F)`) < 0.01 ~ "**",
                                            abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        mutate(across(is.numeric, round, digits=3)) %>% 
        rename("<i>p</i>-value" = "Pr(>F)",
               "R<sup>2</sup>" = "R2",
               "Degree of freedom" = "Df") %>% 
        mutate(`<i>p</i>-value` = format(`<i>p</i>-value`, nsmall = 3)) %>%
        select(c("Degree of freedom", "R<sup>2</sup>", "<i>p</i>-value", " ")) %>% 
        kbl(format = "html", escape = 0) %>%
        kable_styling(full_width = 0, html_font = "serif")

PERMANOVA - Morisita-Horn (stratified)

Distances between samples within each subject. Mean distance between control <-> treatment for each subject

a <- horn_perm_bal %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
                mutate(` ` = case_when(abs(`Pr(>F)`) < 0.001 ~ "***",
                                            abs(`Pr(>F)`) < 0.01 ~ "**",
                                            abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        mutate(across(is.numeric, round, digits=3)) %>% 
        rename("<i>p</i>-value" = "Pr(>F)",
               "R<sup>2</sup>" = "R2",
               "Degree of freedom" = "Df") %>% 
        select(c("R<sup>2</sup>", "<i>p</i>-value", " ")) 

b <- horn_perm_ns %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
                mutate(` ` = case_when(abs(`Pr(>F)`) < 0.001 ~ "***",
                                            abs(`Pr(>F)`) < 0.01 ~ "**",
                                            abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        mutate(across(is.numeric, round, digits=3)) %>% 
        rename("<i>p</i>-value" = "Pr(>F)",
               "R<sup>2</sup>" = "R2",
               "Degree of freedom" = "Df") %>% 
        select(c("R<sup>2</sup>", "<i>p</i>-value", " ")) 

c <- horn_perm_spt %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = case_when(row.names == "lypma" ~ 'lyPMA',
                                     row.names == "benzonase" ~ 'Benzonase',
                                     row.names == "host_zero" ~ 'Host zero',
                                     row.names == "molysis" ~ 'Molysis',
                                     row.names == "qiaamp" ~ 'QIAamp',
                                     row.names == "subject_id" ~ 'Subject id',
                                     row.names == "log10(Final_reads)" ~ 'log10(Final reads)',
                                     row.names == "Residual" ~ 'Residual',
                                     row.names == "Total" ~ 'Total')) %>% column_to_rownames('row.names') %>% 
                mutate(` ` = case_when(abs(`Pr(>F)`) < 0.001 ~ "***",
                                            abs(`Pr(>F)`) < 0.01 ~ "**",
                                            abs(`Pr(>F)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        mutate(across(is.numeric, round, digits=3)) %>% 
        rename("<i>p</i>-value" = "Pr(>F)",
               "R<sup>2</sup>" = "R2",
               "Degree of freedom" = "Df") %>% 
        select(c("R<sup>2</sup>", "<i>p</i>-value", " ")) 


cbind(a, b, c) %>% 
        kbl(format = "html", escape = 0) %>%
        add_header_above(c(" " = 1, "BAL" = 3, "Nasal swab" = 3,"Sputum"= 3)) %>% 
        kable_styling(full_width = 0, html_font = "serif")

Beta diversity boxplot (Horn-Morisita)

Distances between samples within each subject. Mean distance between control <-> treatment for each subject

#distances of betadiversity - boxplots
horn_dist_long <- distance(phyloseq_rel_nz, method="horn") %>% as.matrix() %>% melt_dist() #making long data of distance matrices

#Adding sample type and treatment name. 
#this can be also done by merging metadata into the `horn_dist_long`
names <- data.frame(str_split_fixed(horn_dist_long$iso1, "_", 3))
names2 <- data.frame(str_split_fixed(horn_dist_long$iso2, "_", 3))
horn_dist_long$sample_id_1 <- paste(names$X1, names$X2, sep = "_")
horn_dist_long$method_1 <- ifelse(grepl("lyPMA", horn_dist_long$iso1),"lypma", 
                                         ifelse(grepl("ben", horn_dist_long$iso1),"benzonase", 
                                                ifelse(grepl("host", horn_dist_long$iso1),"host_zero", 
                                                       ifelse(grepl("qia", horn_dist_long$iso1),"qiaamp", 
                                                              ifelse(grepl("moly", horn_dist_long$iso1),"molysis", 
                                                                     "control")))))


#Adding data for iso 2 also should be done
horn_dist_long$sample_id_2 <- paste(names2$X1, names2$X2, sep = "_")
horn_dist_long$method_2 <-ifelse(grepl("lyPMA", horn_dist_long$iso2),"lypma", 
                                        ifelse(grepl("ben", horn_dist_long$iso2),"benzonase", 
                                               ifelse(grepl("host", horn_dist_long$iso2),"host_zero", 
                                                      ifelse(grepl("qia", horn_dist_long$iso2),"qiaamp", 
                                                             ifelse(grepl("moly", horn_dist_long$iso2),"molysis", 
                                                                    "control")))))


#subsetting distances of my interest
horn_dist_long$sample_id_1 <- ifelse(grepl("pos", horn_dist_long$sample_id_1, ignore.case = T),"Mock", 
                                 ifelse(grepl("neg|n_", horn_dist_long$sample_id_1, ignore.case = T),"Neg.",
                                        horn_dist_long$sample_id_1))
horn_dist_long$sample_id_2 <- ifelse(grepl("pos", horn_dist_long$sample_id_2, ignore.case = T),"Mock", 
                                 ifelse(grepl("neg|n_", horn_dist_long$sample_id_2, ignore.case = T),"Neg.",
                                        horn_dist_long$sample_id_2))


path_horn_dist_long_within_sampleid_from_control <- subset(horn_dist_long, horn_dist_long$sample_id_1 == horn_dist_long$sample_id_2) # data within samples

path_horn_dist_long_within_sampleid_from_control <- subset(path_horn_dist_long_within_sampleid_from_control,
                                                           path_horn_dist_long_within_sampleid_from_control$method_1 != path_horn_dist_long_within_sampleid_from_control$method_2) # remove irrelevant association

path_horn_dist_long_within_sampleid_from_control <- subset(path_horn_dist_long_within_sampleid_from_control, (path_horn_dist_long_within_sampleid_from_control$method_1 == "control") + (path_horn_dist_long_within_sampleid_from_control$method_2 == "control") != 0)
path_horn_dist_long_within_sampleid_from_control

path_horn_dist_long_within_sampleid_from_control$treatment <- path_horn_dist_long_within_sampleid_from_control$method_1

path_horn_dist_long_within_sampleid_from_control$treatment <- ifelse(path_horn_dist_long_within_sampleid_from_control$treatment == "control", path_horn_dist_long_within_sampleid_from_control$method_2, path_horn_dist_long_within_sampleid_from_control$treatment) 


#Setting key method
path_horn_dist_long_within_sampleid_from_control$sample_type <- ifelse(grepl("NS", path_horn_dist_long_within_sampleid_from_control$iso1), "Nasal",
                                                                  ifelse(grepl("CFB", path_horn_dist_long_within_sampleid_from_control$iso1), "Sputum",
                                                                         ifelse(grepl("BAL", path_horn_dist_long_within_sampleid_from_control$iso1), "BAL",
                                                                                ifelse(grepl("pos", path_horn_dist_long_within_sampleid_from_control$iso1, ignore.case = T), "Mock",
                                                                                       ifelse(grepl("neg|N_EXT", path_horn_dist_long_within_sampleid_from_control$iso1), "Neg.",NA)))))

#Making a column for baseline (controls, from where?)
path_horn_dist_long_within_sampleid_from_control <- path_horn_dist_long_within_sampleid_from_control %>% 
        mutate(dist_from = case_when(method_1 == "control" ~ iso1,
                                     method_2 == "control" ~ iso2))

dummy <- data.frame(iso1 = path_horn_dist_long_within_sampleid_from_control$dist_from %>% unique,
           iso2 = path_horn_dist_long_within_sampleid_from_control$dist_from %>% unique,
           dist = 0,
           treatment = "Untreated",
           method_1 = "control",
           method_2 = "control"
           )
names <- data.frame(str_split_fixed(dummy$iso1, "_", 3))
names2 <- data.frame(str_split_fixed(dummy$iso2, "_", 3))
dummy$sample_id_1 <- paste(names$X1, names$X2, sep = "_")
#Adding data for iso 2 also should be done
dummy$sample_id_2 <- paste(names2$X1, names2$X2, sep = "_")


#subsetting distances of my interest
dummy$sample_id_1 <- ifelse(grepl("pos", dummy$sample_id_1, ignore.case = T),"Mock", 
                                 ifelse(grepl("neg|n_", dummy$sample_id_1, ignore.case = T),"Neg.",
                                        dummy$sample_id_1))
dummy$sample_id_2 <- ifelse(grepl("pos", dummy$sample_id_2, ignore.case = T),"Mock", 
                                 ifelse(grepl("neg|n_", dummy$sample_id_2, ignore.case = T),"Neg.",
                                        dummy$sample_id_2))
dummy$sample_type <- ifelse(grepl("NS", dummy$iso1), "Nasal",
                            ifelse(grepl("CFB", dummy$iso1), "Sputum",
                                   ifelse(grepl("BAL", dummy$iso1), "BAL",
                                          ifelse(grepl("pos|POS", dummy$iso1, ignore.case = T), "Mock",
                                                 ifelse(grepl("neg|N_EXT", dummy$iso1), "Neg.",NA)))))
dummy <- subset(dummy, !is.na(dummy$sample_type))
path_horn_dist_long_within_sampleid_from_control <- bind_rows(path_horn_dist_long_within_sampleid_from_control, dummy)


path_horn_dist_long_within_sampleid_from_control %>%
        mutate(across(sample_type, factor, levels=c("Mock", "BAL", "Nasal","Sputum")),
               across(treatment, factor, levels=c("Untreated", "lypma", "benzonase", "host_zero","molysis", "qiaamp"))) %>%
ggplot(aes(y = dist, fill = treatment)) +
        geom_boxplot() +
        #scale_fill_manual(values = c(viridis(6)[2:6])) +
        scale_fill_manual(values = c("#e31a1c", "#fb9a99", "#33a02c", "#b2df8a", "#1f78b4", "#a6cee3"), name = "Treatment", labels = c("Untreated", "lyPMA", "Benzonase", "Host zero", "Molysis", "QIAaamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        ylab("Sample pair distances (Horn-Morisita)") +
        theme_classic (base_size = 12, base_family = "serif") + 
        theme(plot.tag = element_text(size = 15),  axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
        facet_wrap(~sample_type, ncol = 5)

LME on H-M distance

Effect size, standard error (SE) and p-value of a statistical test for Bray-Curtis distance from untreated to each treated within subject (BC-distance wihtin subject by treatment~ sample type + treatment + sampetype * treatment + log10 (Final_reads) + (1|subject_id) ).

Hard to interpret, as it does not have proper baseline (control is not there as control does not have distances)

library(lmerTest)
path_horn_dist_long_within_sampleid_from_control$treatment <- factor(path_horn_dist_long_within_sampleid_from_control$treatment,
                                                                     levels = c("Untreated",
                                                                                "lypma",
                                                                                "benzonase",
                                                                                "host_zero",
                                                                                "molysis",
                                                                                "qiaamp"))
path_horn_dist_long_within_sampleid_from_control$sample_type <- factor(path_horn_dist_long_within_sampleid_from_control$sample_type,
                                                                     levels = c("Mock",
                                                                                "BAL",
                                                                                "Nasal",
                                                                                "Sputum"))

tableS9 <- lmer(dist ~ sample_type * treatment + (1|sample_id_1), data = path_horn_dist_long_within_sampleid_from_control %>% subset(., .$sample_type != "Neg.")) %>% summary %>% .$ coefficients %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = gsub("treatment|sample_type", "", row.names)) %>% 
        mutate(row.names = gsub("[:]", " * ", row.names)) %>% 
        mutate(row.names = gsub("host_zero", "Host zero", row.names)) %>% 
        mutate(row.names = gsub("benzonase", "Benzonase", row.names)) %>% 
        mutate(row.names = gsub("lypma", "lyPMA", row.names)) %>% 
        mutate(row.names = gsub("molysis", "Molysis", row.names)) %>% 
        mutate(row.names = gsub("qiaamp", "QIAamp", row.names)) %>% 
        mutate(across(is.numeric, round, digits=3)) %>% 
        column_to_rownames("row.names") %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.001 ~ "***",
                                            abs(`Pr(>|t|)`) < 0.01 ~ "**",
                                            abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rename("<i>p</i>-value" = "Pr(>|t|)",
               SE = "Std. Error") %>%
        mutate("Effect size (95% CI)" = paste(round(Estimate, 2) %>% format(nsmall = 2), 
                                " (", 
                                round(Estimate - 1.96 * SE, 2) %>% format(nsmall = 2),
                                ", ",
                                round(Estimate + 1.96 * SE, 2) %>% format(nsmall = 2),
                                ")", 
                                sep = ""),
               "<i>p</i>-value" = round(`<i>p</i>-value`, 3)) %>%
        select(c("Effect size (95% CI)", "<i>p</i>-value", " ")) %>% 
        kbl(format = "html", escape = 0) %>%
        kable_styling(full_width = 0, html_font = "serif")


tableS9

LME on H-M distance (stratified)

path_horn_dist_long_within_sampleid_from_control$treatment <- factor(path_horn_dist_long_within_sampleid_from_control$treatment,
                                                                     levels = c("Untreated",
                                                                                "lypma",
                                                                                "benzonase",
                                                                                "host_zero",
                                                                                "molysis",
                                                                                "qiaamp"))
path_horn_dist_long_within_sampleid_from_control$sample_type <- factor(path_horn_dist_long_within_sampleid_from_control$sample_type,
                                                                     levels = c("BAL",
                                                                                "Nasal",
                                                                                "Sputum"))

a <- lm(dist ~ treatment, data = path_horn_dist_long_within_sampleid_from_control %>% subset(., .$sample_type == "BAL")) %>% summary %>% .$ coefficients %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = gsub("host_zero", "Host zero", row.names)) %>% 
        mutate(row.names = gsub("benzonase", "Benzonase", row.names)) %>% 
        mutate(row.names = gsub("lypma", "lyPMA", row.names)) %>% 
        mutate(row.names = gsub("molysis", "Molysis", row.names)) %>% 
        mutate(row.names = gsub("qiaamp", "QIAamp", row.names)) %>% 
        mutate(across(is.numeric, round, digits=3)) %>% 
        mutate(`row.names` = gsub("treatment|sample_type", "", `row.names`)) %>% 
        column_to_rownames("row.names") %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.001 ~ "***",
                                            abs(`Pr(>|t|)`) < 0.01 ~ "**",
                                            abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rename("<i>p</i>-value" = "Pr(>|t|)",
               SE = "Std. Error") %>%
        mutate("Effect size (95% CI)" = paste(round(Estimate, 2) %>% format(nsmall = 2), 
                                " (", 
                                round(Estimate - 1.96 * SE, 2) %>% format(nsmall = 2),
                                ", ",
                                round(Estimate + 1.96 * SE, 2) %>% format(nsmall = 2),
                                ")", 
                                sep = ""),
               "<i>p</i>-value" = round(`<i>p</i>-value`, 3)) %>%
        select(c("Effect size (95% CI)", "<i>p</i>-value", " ")) 

b <- lm(dist ~ treatment, data = path_horn_dist_long_within_sampleid_from_control %>% subset(., .$sample_type == "Nasal")) %>% summary %>% .$ coefficients %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = gsub("host_zero", "Host zero", row.names)) %>% 
        mutate(row.names = gsub("benzonase", "Benzonase", row.names)) %>% 
        mutate(row.names = gsub("lypma", "lyPMA", row.names)) %>% 
        mutate(row.names = gsub("molysis", "Molysis", row.names)) %>% 
        mutate(row.names = gsub("qiaamp", "QIAamp", row.names)) %>% 
        mutate(across(is.numeric, round, digits=3)) %>% 
        mutate(`row.names` = gsub("treatment|sample_type", "", `row.names`)) %>% 
        column_to_rownames("row.names") %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.001 ~ "***",
                                            abs(`Pr(>|t|)`) < 0.01 ~ "**",
                                            abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rename("<i>p</i>-value" = "Pr(>|t|)",
               SE = "Std. Error") %>%
        mutate("Effect size (95% CI)" = paste(round(Estimate, 2) %>% format(nsmall = 2), 
                                " (", 
                                round(Estimate - 1.96 * SE, 2) %>% format(nsmall = 2),
                                ", ",
                                round(Estimate + 1.96 * SE, 2) %>% format(nsmall = 2),
                                ")", 
                                sep = ""),
               "<i>p</i>-value" = round(`<i>p</i>-value`, 3)) %>%
        select(c("Effect size (95% CI)", "<i>p</i>-value", " ")) 

c <- lm(dist ~ treatment, data = path_horn_dist_long_within_sampleid_from_control %>% subset(., .$sample_type == "Sputum")) %>% summary %>% .$ coefficients %>% data.frame(check.names = F) %>% rownames_to_column('row.names') %>% 
        mutate(row.names = gsub("host_zero", "Host zero", row.names)) %>% 
        mutate(row.names = gsub("benzonase", "Benzonase", row.names)) %>% 
        mutate(row.names = gsub("lypma", "lyPMA", row.names)) %>% 
        mutate(row.names = gsub("molysis", "Molysis", row.names)) %>% 
        mutate(row.names = gsub("qiaamp", "QIAamp", row.names)) %>% 
        mutate(across(is.numeric, round, digits=3)) %>% 
        mutate(`row.names` = gsub("treatment|sample_type", "", `row.names`)) %>% 
        column_to_rownames("row.names") %>% 
        mutate(` ` = case_when(abs(`Pr(>|t|)`) < 0.001 ~ "***",
                                            abs(`Pr(>|t|)`) < 0.01 ~ "**",
                                            abs(`Pr(>|t|)`) < 0.05 ~ "*",
                               .default = " ")) %>% 
        rename("<i>p</i>-value" = "Pr(>|t|)",
               SE = "Std. Error") %>%
        mutate("Effect size (95% CI)" = paste(round(Estimate, 2) %>% format(nsmall = 2), 
                                " (", 
                                round(Estimate - 1.96 * SE, 2) %>% format(nsmall = 2),
                                ", ",
                                round(Estimate + 1.96 * SE, 2) %>% format(nsmall = 2),
                                ")", 
                                sep = ""),
               "<i>p</i>-value" = round(`<i>p</i>-value`, 3)) %>%
        select(c("Effect size (95% CI)", "<i>p</i>-value", " "))

cbind(a, b, c) %>% 
        kbl(format = "html", escape = 0) %>%
        add_header_above(c(" " = 1, "BAL" = 3, "Nasal swab" = 3,"Sputum"= 3)) %>% 
        kable_styling(full_width = 0, html_font = "serif")

Factors affecting DA taxa (main text)

filt_maaslin_all <- read.csv("data/filt_maaslin_all.csv")
filt_maaslin_interaction <- read.csv("data/filt_maaslin_interaction.csv")
filt_fit_data_bal <- read.csv("data/filt_fit_data_bal.csv")
filt_fit_data_spt <- read.csv("data/filt_fit_data_spt.csv")
filt_fit_data_ns <- read.csv("data/filt_fit_data_ns.csv")
filt_fit_data_pos <- read.csv("data/filt_fit_data_pos.csv")

cat("Factors affecting DA taxa (q<0.1)")

## Factors affecting DA taxa (q<0.1)

filt_maaslin_all %>% subset(., .$qval < 0.1 ) %>% .$metadata %>% table

## .
##         benzonase         host_zero log10.Final_reads             lypma 
##                34                19                85                31 
##           molysis            qiaamp       sample_type 
##                31                22               118

Volcano plot

Volcano plot of sequencing depth adjusted differential abundance of taxa by each treatment

#Making significance table for figure
        # Define a function to make species names italicized
# Make a significance table for each figure (top 20 taxa)
species_italic <- function(data) {
  names <- gsub("_", " ", rownames(data))
  names <- gsub("[]]|[[]", "", names)
  names <- gsub(" sp", " sp.", names)
  names <- gsub(" sp.", "* sp.", names)
  names <- gsub(" group", "* group.", names)
  names <- ifelse(grepl("[*]", names), paste("*", names, sep = ""), paste("*", names, "*", sep = ""))
  rownames(data) <- names
  data
}

make_sig_table <- function(data) {
  sig_data <- spread(data[order(data$qval), c("feature", "metadata", "qval")], metadata, qval)
  sig_data$min <- apply(sig_data, 1, FUN = min)
  sig_data <- sig_data[order(sig_data$min),] %>% select("feature", "lypma", "benzonase", "host_zero", "molysis", "qiaamp") %>% .[1:20,]
  sig_data[["feature"]] <- ifelse(sig_data[["feature"]] == "X.Collinsella._massiliensis", "[Collinsella]_massiliensis", sig_data[["feature"]])
  sig_data_italic <- sig_data %>% rownames_to_column(var = "-") %>%
          column_to_rownames(var = "feature") %>% species_italic %>% select(-c("-")) %>%
          rename(lyPMA = lypma,  Benzonase = benzonase, `Host zero` = host_zero, Molysis = molysis, QIAamp = qiaamp)
  sig_data_sig <- ifelse(sig_data_italic < 0.1, "*", NA) %>% data.frame(check.names = F)
  return(list(data = sig_data, data_italic = sig_data_italic, data_sig = sig_data_sig))
}
filt_fit_data_pos <- make_sig_table(filt_fit_data_pos)
filt_fit_data_bal <- make_sig_table(filt_fit_data_bal)
filt_fit_data_ns <- make_sig_table(filt_fit_data_ns)
filt_fit_data_spt <- make_sig_table(filt_fit_data_spt)

filt_pos_sig <- subset_taxa(subset_samples(phyloseq_unfiltered$phyloseq_rel, sample_type == "Mock"),
                                       taxa_names(subset_samples(phyloseq_unfiltered$phyloseq_rel, sample_type == "Mock")) %in% filt_fit_data_pos$data$feature)

filt_fit_data_pos$rel <- cbind(filt_pos_sig %>% otu_table %>% t, filt_pos_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>%  data.frame(check.names = F) %>% 
        .[row.names(filt_fit_data_pos$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")

filt_spt_sig <- subset_taxa(subset_samples(phyloseq_unfiltered$phyloseq_rel, sample_type == "Sputum"), taxa_names(subset_samples(phyloseq_unfiltered$phyloseq_rel, sample_type == "Sputum")) %in% filt_fit_data_spt$data$feature)

filt_fit_data_spt$rel <- cbind(filt_spt_sig %>% otu_table %>% t, filt_spt_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>%  data.frame(check.names = F) %>% 
        .[row.names(filt_fit_data_spt$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")

filt_ns_sig <- subset_taxa(subset_samples(phyloseq_unfiltered$phyloseq_rel, sample_type == "Nasal"),
                                       taxa_names(subset_samples(phyloseq_unfiltered$phyloseq_rel, sample_type == "Nasal")) %in% filt_fit_data_ns$data$feature)

filt_fit_data_ns$rel <- cbind(filt_ns_sig %>% otu_table %>% t, filt_ns_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>%  data.frame(check.names = F) %>% 
        .[row.names(filt_fit_data_ns$data_italic),] %>%  mutate_all(~na_if(., 0)) %>% rownames_to_column("feature")

filt_bal_sig <- subset_taxa(subset_samples(phyloseq_unfiltered$phyloseq_rel, sample_type == "BAL"),
                                       taxa_names(subset_samples(phyloseq_unfiltered$phyloseq_rel, sample_type == "BAL")) %in% filt_fit_data_bal$data$feature)

filt_fit_data_bal$rel <- cbind(filt_bal_sig %>% otu_table %>% t, filt_bal_sig %>% sample_data) %>% group_by(treatment) %>% summarise_if(is.numeric, mean, na.rm = TRUE) %>% .[, 1:21] %>% column_to_rownames(., "treatment") %>% t () %>% species_italic() %>%  data.frame(check.names = F) %>%
        .[row.names(filt_fit_data_bal$data_italic),] %>%
        mutate_all(~na_if(., 0)) %>% rownames_to_column("feature") %>% subset(., !grepl("NA", .$feature))

#Volcano plot

ggplot(filt_maaslin_all, aes(y = -log10(qval), x = coef, col = metadata)) +
        theme_classic(base_family = "serif") +
        #labs(tag = "A") +
        geom_point(size = 2) +
        xlab("MaAslin coefficient") +
        ylab("-log<sub>10</sub>(*q*-value)") +
        ylim(c(-1, 35)) +
        geom_hline(yintercept = 1, col = "gray") +
        geom_vline(xintercept = 0, col = "gray") +
        annotate(family = "serif",
                 geom='richtext',
                 x=0, y=80,
                 label = "<i>q</i>-value = 0.1, fold-change = 0") +
        theme(legend.position = "top", axis.title.y = ggtext::element_markdown(), legend.text = element_markdown()) +
        scale_color_manual(values = c("#4daf4a",  "#984ea3", "#f781bf", "#377eb8", "#ff7f00", "#ffff33", "#a65628"),
                           breaks = c("log10.Final_reads", "sample_type", "lypma", "benzonase", "host_zero",  "molysis", "qiaamp"), 
                           labels = c("log<sub>10</sub>(Final reads)", "Sample type", "lyPMA", "Benzonase", "Host zero",  "Molysis", "QIAamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        guides(col = guide_legend(title = "Covariates", title.position = "top", nrow = 1))

Balloon plot

Mean relative abundance of top 20 significant taxa identified by differential abundance analysis by sample type. (A) Mock community, (B) bronchoalveolor lavage, (C) nasal swabs, and (D) sputum.

f5a <- merge(filt_fit_data_pos$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      filt_fit_data_pos$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(filt_fit_data_pos$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
        mutate(sig = case_when(sig < 0.1 ~ "< 0.1",
                               .default = "> 0.1"),
               feature = case_when(feature =="*Saccharomyces cerevisiae x Saccharomyces kudriavzevii*" ~ "*S. cerevisiae* X *S. kudriavzevii*",
                               .default = feature)) %>%
#Baloon plot
        ggballoonplot(size = "value", y = "feature", x= "treatment", fill = "sig") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "A") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown(size = 7))  +
        #Adding significance asterisks
        scale_fill_manual(values = c("red", "grey"), aes(y = feature,
                      x = treatment,
                      label = sig)) +
        guides(fill = guide_legend(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top"),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 1),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         )
#ffff33 qia

f5b <- merge(filt_fit_data_bal$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      filt_fit_data_bal$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(filt_fit_data_bal$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
        mutate(sig = case_when(is.na(qval) ~ "> 0.1",
                               qval < 0.1 ~ "< 0.1",
                               .default = "> 0.1")) %>%

#Baloon plot
        ggballoonplot(size = "value", y = "feature", x= "treatment", fill = "sig") +
        theme_classic(base_family = "serif") +
        #colors for qvalues
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "B") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown(size = 7))  +
        
        scale_fill_manual(values = c("grey", "red"), aes(y = feature,
                      x = treatment,
                      label = sig)) +
        guides(fill = guide_legend(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top"),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 1),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         )

f5c <- merge(filt_fit_data_ns$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      filt_fit_data_ns$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(filt_fit_data_ns$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
        mutate(sig = case_when(sig < 0.1 ~ "< 0.1",
                               .default = "> 0.1")) %>%
#Baloon plot
        ggballoonplot(size = "value", y = "feature", x= "treatment", fill = "sig") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "C") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown(size = 7))  +
        #Adding significance asterisks
        scale_fill_manual(values = c("red", "grey"), aes(y = feature,
                      x = treatment,
                      label = sig)) +
        guides(fill = guide_legend(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top"),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 1),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         )

f5d <- merge(filt_fit_data_spt$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      filt_fit_data_spt$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        merge(filt_fit_data_spt$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
        mutate(sig = case_when(sig < 0.1 ~ "< 0.1",
                               .default = "> 0.1")) %>%
#Baloon plot
        ggballoonplot(size = "value", y = "feature", x= "treatment", fill = "sig") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Species") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown(size = 7),
              axis.title.x = element_text(margin = margin(t = 20)))  +
        #Adding significance asterisks
        scale_fill_manual(values = c("red", "grey"), aes(y = feature,
                      x = treatment,
                      label = sig)) +
        guides(fill = guide_legend(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top"),
               size = guide_legend(title = "Relative abundance",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 1),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         )
ggarrange(f5a, f5b, f5c, f5d, align = "hv", common.legend = T)

MaAslin function - volcano plot (Fig. S9)

Figure S9. Volcano plot of sequencing depth adjusted differential abundance of function by each treatment

#Volcano plot

figS9 <- ggplot(f_maaslin_all, aes(y = -log10(qval), x = coef, col = metadata)) +
        theme_classic(base_family = "serif") +
        #labs(tag = "A") +
        geom_point(size = 2) +
        xlab("MaAslin coefficient") +
        ylab("-log<sub>10</sub>(*q*-value)") +
        geom_hline(yintercept = 1, col = "gray") +
        geom_vline(xintercept = 0, col = "gray") +
        annotate(family = "serif",
                 geom='richtext',
                 x=0, y=80,
                 label = "<i>q</i>-value = 0.1, fold-change = 0") +
        theme(legend.position = "top", axis.title.y = ggtext::element_markdown()) +
        scale_color_manual(values = c("#4daf4a",  "#984ea3", "#f781bf", "#377eb8", "#ff7f00", "#ffff33", "#a65628"),
                           breaks = c("log10.Final_reads", "sample_type", "lypma", "benzonase", "host_zero",  "molysis", "qiaamp"), 
                           labels = c("log10(Final reads)", "Sample type", "lyPMA", "Benzonase", "Host zero",  "Molysis", "QIAamp")) + #color using https://colorbrewer2.org/#type=qualitative&scheme=Set1&n=6
        guides(col = guide_legend(title = "Covariates", title.position = "top", nrow = 2))


png(file = "/Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/7_Manuscripts/2022_MGK_Host_Depletion/Figures/FigureS9.png",   # The directory you want to save the file in
    width = 180, # The width of the plot in inches
    height = 90, # The height of the plot in inches
    units = "mm",
    res = 600
) #fixing multiple page issue

figS9
# alpha diversity plots
#ggarrange(f4ad, ggarrange(f4e, f4f, ncol = 2),
#          ncol = 1) # alpha diversity plots

dev.off()

## quartz_off_screen 
##                 2

figS9

Function baloon plot (Fig. S10)

Figure S10. Mean relative abundance of top 20 significant function identified by differential abundance analysis using MaAsLin. Analyses were stratified by sample type. (A) Mock community, (B) bronchoalveolor lavage, (C) nasal swabs, and (D) sputum. Statistical significances were noted at the level of q-value < 0.1

f5a <- merge(f_fit_data_pos$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      f_fit_data_pos$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(f_fit_data_pos$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
        mutate(sig = case_when(sig < 0.1 ~ "< 0.1",
                               .default = "> 0.1"),
               value = value * 1000000) %>%
#Baloon plot
        ggballoonplot(size = "value", y = "feature", x= "treatment", fill = "sig") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Function") +
        labs(tag = "A") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
                axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown(size = 7),
              axis.title.x = element_text(margin = margin(t = 20)))  +
        #Adding significance asterisks
        scale_fill_manual(values = c("red", "grey"), aes(y = feature,
                      x = treatment,
                      label = sig)) +
        guides(fill = guide_legend(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top"),
               size = guide_legend(title = "Copies per million",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 1),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         )
#ffff33 qia

f5b <- merge(f_fit_data_bal$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      f_fit_data_bal$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(f_fit_data_bal$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
        mutate(sig = case_when(sig < 0.1 ~ "< 0.1",
                               .default = "> 0.1"),
               value = value * 1000000) %>%

#Baloon plot
        ggballoonplot(size = "value", y = "feature", x= "treatment", fill = "sig") +
        theme_classic(base_family = "serif") +
        #colors for qvalues
        xlab("Experimental group") +
        ylab("Function") +
        labs(tag = "B") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown(size = 7),
              axis.title.x = element_text(margin = margin(t = 20)))  +
        
        scale_fill_manual(values = c("grey", "red"), aes(y = feature,
                      x = treatment,
                      label = sig)) +
        guides(fill = guide_legend(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top"),
               size = guide_legend(title = "Copies per million",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 1),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         )

f5c <- merge(f_fit_data_ns$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      f_fit_data_ns$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        
        merge(f_fit_data_ns$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
        mutate(sig = case_when(sig < 0.1 ~ "< 0.1",
                               .default = "> 0.1"),
               value = value * 1000000) %>%
#Baloon plot
        ggballoonplot(size = "value", y = "feature", x= "treatment", fill = "sig") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Function") +
        labs(tag = "C") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown(size = 7),
              axis.title.x = element_text(margin = margin(t = 20)))  +
        #Adding significance asterisks
        scale_fill_manual(values = c("red", "grey"), aes(y = feature,
                      x = treatment,
                      label = sig)) +
        guides(fill = guide_legend(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top"),
               size = guide_legend(title = "Copies per million",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 1),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         )

f5d <- merge(f_fit_data_spt$rel %>%
              gather(treatment,
                     value,
                     Untreated:QIAamp,
                     factor_key=TRUE),
      f_fit_data_spt$data_italic %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     qval,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
      by.x = c('feature', 'treatment'),
      by.y = c('feature', 'treatment'),
      all = T) %>%
        merge(f_fit_data_spt$data_sig %>%
              rownames_to_column("feature") %>%
              gather(treatment,
                     sig,
                     lyPMA:QIAamp,
                     factor_key=TRUE),
              by.x = c('feature', 'treatment'),
              by.y = c('feature', 'treatment'),
              all = T) %>%
        mutate(sig = case_when(sig < 0.1 ~ "< 0.1",
                               .default = "> 0.1"),
               value = value * 1000000) %>%
#Baloon plot
        ggballoonplot(size = "value", y = "feature", x= "treatment", fill = "sig") +
        
        theme_classic(base_family = "serif") +
        #colors for qvalues
        gradient_fill(c("#006d2c", "#edf8fb")) +
        xlab("Experimental group") +
        ylab("Function") +
        labs(tag = "D") +
        theme(panel.grid.major = element_line(colour = "grey"),
              legend.position = "top",
              axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
              #Element markdown for taxa name italicizing
              axis.text.y = ggtext::element_markdown(size = 7),
              axis.title.x = element_text(margin = margin(t = 20)))  +
        #Adding significance asterisks
        scale_fill_manual(values = c("red", "grey"), aes(y = feature,
                      x = treatment,
                      label = sig)) +
        guides(fill = guide_legend(title = c(expression(paste(italic("q"),
                                                       "-value",
                                                       sep = ""))),
                                   title.position = "top"),
               size = guide_legend(title = "Copies per millio",
                                   title.position = "top",
                                   order = 1,
                                   nrow = 1),
               ) + 
        scale_x_discrete(labels=c("control" = "Untreated",
                                  "lypma" = "lyPMA",
                                  "benzonase" = "Benzonase",
                                  "host_zero" = "Host-zero",
                                  "molysis" = "Molysis",
                                  "qiaamp" = "QIAamp")
                         )

figS10 <- ggarrange(f5a, f5b, f5c, f5d, align = "hv", common.legend = T)

figS10

png(file = "/Users/minsikkim/Dropbox (Partners HealthCare)/Project_SICAS2_microbiome/7_Manuscripts/2022_MGK_Host_Depletion/Figures/FigureS10.png",   # The directory you want to save the file in
    width = 180, # The width of the plot in inches
    height = 170, # The height of the plot in inches
    units = "mm",
    res = 600
) #fixing multiple page issue

figS10
# alpha diversity plots
#ggarrange(f4ad, ggarrange(f4e, f4f, ncol = 2),
#          ncol = 1) # alpha diversity plots

dev.off()

## quartz_off_screen 
##                 2

Rarefaction stratified (BAL)

args=(commandArgs(TRUE))

if(length(args)==0){
    print("No arguments supplied")
    
    dist_type <- 5      # 1:"horn"; 5: "jaccard"; 4: "unwt-unifrac"
    data_full <- 1      # 1: full sample; 0: a subsample of 50 subjects
    i_seed <- seed         # seed for generating different subsamples
    
} else{
    for(i in 1:length(args)){
        eval(parse(text=args[[i]]))
    }
}

## [1] "No arguments supplied"

if (dist_type==5) {
    dist_method = "jaccard"
} else if (dist_type==4) {
    dist_method = "unwt-unifrac"
} else if (dist_type==1) {
    dist_method = "horn"
}  


###############################################
# read in data
###############################################

# read in otu table, taxonomy info
#dat <- read_biom("rawData/otu_table.biom")
#otu_table <- as.data.frame(as.matrix(biom_data(dat)))
#taxonomy <- observation_metadata(dat)
#crohn.tree <- read.tree("rawData/insertion_tree.relabelled.tre")

otu_table_col <- as.data.frame(t(otu_table(phyloseq$phyloseq_count %>%
                                                   subset_samples(sample_type %in% c("BAL"))))) #make samples rows, OTUs columns
otu_table_col <- otu_table_col[,colSums(otu_table_col>0)>0] # Removing data with zeros

# read in metadata
metadata <- sample_data(phyloseq$phyloseq_count %>%
                                                   subset_samples(sample_type %in% c("BAL"))) %>% data.frame()

# create dichotomous variable for diagnosis
#metadata <- cbind(metadata,ifelse(metadata$diagnosis=="CD" | metadata$diagnosis=="IC" | metadata$diagnosis=="UC", 1, 0))
#colnames(metadata)[61] <- "disease" 

# sort metadata to match otu table
meta.sorted <- metadata[order(match(metadata$baylor_other_id, rownames(otu_table_col))),]

# remove duplicates at the end
#meta.site.dup <- subset(meta.sorted, collection=="RISK")
#meta.sub.dup <- subset(meta.site.dup, type_sample=="biopsy" & biopsy_location=="Rectum")

# want subset with only biopsy samples, no follow-up measurements, no missing values in key covariates
#no.follow <- meta.sorted[!duplicated(meta.sorted$anonymized_name),]
#meta.site <- subset(no.follow, collection=="RISK")
#meta.sub <- subset(meta.site, type_sample=="biopsy" & biopsy_location=="Rectum")

# match otu table with cleaned meta data
otu.sub <- subset(otu_table_col, rownames(otu_table_col) %in% metadata$baylor_other_id) %>% filter(.,)

###Prevalence filtration & decontam was not processed as phyloseq was alread filtered
# filter otu table 

pa.table <- 1*(otu.sub > 0) #Presence absence table
#otu.sub.f <- otu.sub[, -which(colSums(pa.table)<5)]  #Prevalence filtering
#dim(otu.sub)
#dim(otu.sub.f)    # 16 taxa were removed

# remove obs with < 5000 reads

otu.full <- otu.sub[rowSums(otu.sub)>5000,] #samples with too low reads were removed
#Subsetting samples with effective reads
meta.sub <- subset(meta.sorted, meta.sorted$baylor_other_id %in% rownames(otu.full))
dim(otu.sub)

## [1] 30 88

dim(otu.full) #14 samples were thrown out

## [1] 27 88

min(rowSums(otu.full))

## [1] 10612

meta.full <- subset(meta.sub, meta.sub$baylor_other_id %in% rownames(otu.full))
#Making data of t-tests - Not usefull for my analysis
#Anyway ran for differences between control vs treatment
p.t.full = t.test(rowSums(otu.full) ~ meta.full$control, var.equal=FALSE)$p.value
p.wilcox.full = wilcox.test(rowSums(otu.full) ~ meta.full$control)$p.value

#1st. We will check all the samples, non-rarefied

########################
# prepare data for anal
########################

if (data_full==1) {
    otu.tab = otu.full
    meta.data = meta.full
} else {
    otu.tab = otu.sam10
    meta.data = meta.sam10
}
dim(otu.tab)

## [1] 27 88

dim(meta.data)

## [1] 27 66

########################
# permanova-F
########################
#LDM-F: averaging the F statistics over all rarefactions, 


if (data_full==1) {
    formula = otu.tab ~ treatment
} else {
    formula = otu.tab ~ treatment
}

res.F <- permanovaFL(formula = formula,
                     data=meta.sub,dist.method = dist_method,
                     scale.otu.table = FALSE, n.rarefy=100, n.rej.stop=100, seed=seed, perm.between.type = "none",
                     perm.within.type = "free", cluster.id = "subject_id")

## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## PERMANOVA stopped at permutation 1100

        #Permutation within requires the same sie of clusters - as some of samples have failed sequencing this option is not valid.
# subject_id was added as a fixed term

p_permanovaF <- res.F$p.permanova


# -----------------------------
# dist for permanova-D2 (1-100)
# -----------------------------

n_rarefy = 100
label = colnames(otu.full)

n_sam_new = nrow(otu.tab)
dist_sq_all = array(0, dim=c(n_sam_new, n_sam_new, n_rarefy))

set.seed(seed)

for (r in 1:n_rarefy) {
    otu_table_rarefy <- Rarefy(otu.tab)$otu.tab.rff
    otu_table_rarefy1 <- (otu_table_rarefy>0)*1
    colnames(otu_table_rarefy1) = label
    colnames(otu_table_rarefy) = label
    if (dist_type == 5 ) {
        dist <- vegdist(x=otu_table_rarefy1, method=dist_method)
    } else if (dist_type == 4) {
        dist <- GUniFrac(otu_table_rarefy1, crohn.tree, alpha=c(1))$unifrac[,,"d_UW"]
    }
    dist <- as.matrix(dist)
    
    if (r==1) {
        dist_sq_all[,,1] = dist^2
    } else {
        dist_sq_all[,,r] = dist_sq_all[,,r-1] + dist^2
    }
}

# --------------------------
# dist for permanova-D2-A
# --------------------------
#LDM-A: averaging the RSS terms over all rarefactions, 
if (dist_type == 5 ) {
    res <- jaccard.mean(otu.tab)
    dist_sq_inf_o1 = res$jac.mean.sq.o1
    dist_sq_inf_o2 = res$jac.mean.sq.o2
} else if (dist_type == 4) {
    res <- unifrac.mean(otu.tab, crohn.tree)
    dist_sq_inf_o1 = res$unifrac.mean.sq.o1
    dist_sq_inf_o2 = res$unifrac.mean.sq.o2
}
#Negative controls were having issues - 0 reads after rarefying. These OTU table with zeros caused problem in data processing, therefore the Neg.s were removed from rarefaction analysis


#########################
# permanova-D2, D2-A
#########################

if (data_full==1) {
        
        formula = dist_sq_r ~ treatment

} else {
    formula = dist_sq_r ~ treatment
}

p_permanovaD2 = NULL
p_permanovaD2A2 = NULL
p_permanovaD2A1 = NULL

for (r in c(1:10, 50, 100, 100^2, 100^3)) { # "100^2" corresponds to D2-A2, "100^3" corresponds to D2-A1
    cat("number of rarefaction:", r, "\n")

    if (r == 100^2) {
        dist_sq_r = dist_sq_inf_o2
    } else if (r == 100^3) {
        dist_sq_r = dist_sq_inf_o1
    } else {
        dist_sq_r = dist_sq_all[,,r]/r
    }
    res_permanovaD2 = permanovaFL(formula = formula, 
                                  data = meta.data, n.perm.max=5000, 
                                  square.dist=FALSE, n.rarefy=0, n.rej.stop=100, seed=seed,
                                  )
    
    if (r == 100^2) {
        p_permanovaD2A2 = res_permanovaD2$p.permanova
    } else if (r == 100^3) {
        p_permanovaD2A1 = res_permanovaD2$p.permanova
    } else {
        p_permanovaD2 = c(p_permanovaD2, res_permanovaD2$p.permanova)
    }
}

## number of rarefaction: 1 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 2 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 3 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 4 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 5 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 6 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 7 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 8 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 9 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 10 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 50 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 100 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 10000 
## permutations: 100 
## PERMANOVA stopped at permutation 100 
## number of rarefaction: 1e+06 
## permutations: 100 
## PERMANOVA stopped at permutation 100

#tab = c(p_permanovaF, p_permanovaD2, p_permanovaD2A2, p_permanovaD2A1, i_seed)
#write.table(t(tab), file = paste("", data_full, "_dist", dist_type, ".txt", sep=""), 
#            quote=FALSE, append=TRUE, row.names=FALSE, col.names=FALSE)


p_permanovaF

##    rarefy1    rarefy2    rarefy3    rarefy4    rarefy5    rarefy6    rarefy7 
## 0.10090909 0.09909091 0.09727273 0.09909091 0.10090909 0.10000000 0.09818182 
##    rarefy8    rarefy9   rarefy10   rarefy11   rarefy12   rarefy13   rarefy14 
## 0.09636364 0.09636364 0.09636364 0.09636364 0.09636364 0.09636364 0.09727273 
##   rarefy15   rarefy16   rarefy17   rarefy18   rarefy19   rarefy20   rarefy21 
## 0.09727273 0.10000000 0.09909091 0.09909091 0.09909091 0.10090909 0.10000000 
##   rarefy22   rarefy23   rarefy24   rarefy25   rarefy26   rarefy27   rarefy28 
## 0.10000000 0.10000000 0.09909091 0.10000000 0.10000000 0.10000000 0.10000000 
##   rarefy29   rarefy30   rarefy31   rarefy32   rarefy33   rarefy34   rarefy35 
## 0.09909091 0.09909091 0.09909091 0.09909091 0.09909091 0.09909091 0.09909091 
##   rarefy36   rarefy37   rarefy38   rarefy39   rarefy40   rarefy41   rarefy42 
## 0.09909091 0.09909091 0.09909091 0.09909091 0.09909091 0.10000000 0.10000000 
##   rarefy43   rarefy44   rarefy45   rarefy46   rarefy47   rarefy48   rarefy49 
## 0.10000000 0.10000000 0.10000000 0.09909091 0.09909091 0.09909091 0.09909091 
##   rarefy50   rarefy51   rarefy52   rarefy53   rarefy54   rarefy55   rarefy56 
## 0.09818182 0.09818182 0.09818182 0.09818182 0.09818182 0.09818182 0.09818182 
##   rarefy57   rarefy58   rarefy59   rarefy60   rarefy61   rarefy62   rarefy63 
## 0.09818182 0.09818182 0.09818182 0.09818182 0.09818182 0.09818182 0.09909091 
##   rarefy64   rarefy65   rarefy66   rarefy67   rarefy68   rarefy69   rarefy70 
## 0.09909091 0.09909091 0.09909091 0.09818182 0.09818182 0.09909091 0.09909091 
##   rarefy71   rarefy72   rarefy73   rarefy74   rarefy75   rarefy76   rarefy77 
## 0.09909091 0.09909091 0.09909091 0.09818182 0.09818182 0.09818182 0.09909091 
##   rarefy78   rarefy79   rarefy80   rarefy81   rarefy82   rarefy83   rarefy84 
## 0.09909091 0.09818182 0.09818182 0.09909091 0.09818182 0.09818182 0.09909091 
##   rarefy85   rarefy86   rarefy87   rarefy88   rarefy89   rarefy90   rarefy91 
## 0.10000000 0.10000000 0.10000000 0.10000000 0.10000000 0.10000000 0.10000000 
##   rarefy92   rarefy93   rarefy94   rarefy95   rarefy96   rarefy97   rarefy98 
## 0.10000000 0.10000000 0.10000000 0.10000000 0.10000000 0.10000000 0.10000000 
##   rarefy99  rarefy100 
## 0.10000000 0.10000000

p_permanovaD2

##  [1] 1 1 1 1 1 1 1 1 1 1 1 1

p_permanovaD2A2

## [1] 1

p_permanovaD2A1

## [1] 1

res_permanovaD2$R.squared

## [1] 0.08734907

Rarefaction stratified (Nasal swab)

As a sanity check, rarefaction based analysis was conducted using LDM

if(length(args)==0){
    print("No arguments supplied")
    
    dist_type <- 5      # 1:"horn"; 5: "jaccard"; 4: "unwt-unifrac"
    data_full <- 1      # 1: full sample; 0: a subsample of 50 subjects
    i_seed <- seed         # seed for generating different subsamples
    
} else{
    for(i in 1:length(args)){
        eval(parse(text=args[[i]]))
    }
}

## [1] "No arguments supplied"

if (dist_type==5) {
    dist_method = "jaccard"
} else if (dist_type==4) {
    dist_method = "unwt-unifrac"
} else if (dist_type==1) {
    dist_method = "horn"
}  


###############################################
# read in data
###############################################

# read in otu table, taxonomy info
#dat <- read_biom("rawData/otu_table.biom")
#otu_table <- as.data.frame(as.matrix(biom_data(dat)))
#taxonomy <- observation_metadata(dat)
#crohn.tree <- read.tree("rawData/insertion_tree.relabelled.tre")

otu_table_col <- as.data.frame(t(otu_table(phyloseq$phyloseq_count %>%
                                                   subset_samples(sample_type %in% c("Nasal"))))) #make samples rows, OTUs columns
otu_table_col <- otu_table_col[,colSums(otu_table_col>0)>0] # Removing data with zeros

# read in metadata
metadata <- sample_data(phyloseq$phyloseq_count %>%
                                                   subset_samples(sample_type %in% c("Nasal"))) %>% data.frame()

# create dichotomous variable for diagnosis
#metadata <- cbind(metadata,ifelse(metadata$diagnosis=="CD" | metadata$diagnosis=="IC" | metadata$diagnosis=="UC", 1, 0))
#colnames(metadata)[61] <- "disease" 

# sort metadata to match otu table
meta.sorted <- metadata[order(match(metadata$baylor_other_id, rownames(otu_table_col))),]

# remove duplicates at the end
#meta.site.dup <- subset(meta.sorted, collection=="RISK")
#meta.sub.dup <- subset(meta.site.dup, type_sample=="biopsy" & biopsy_location=="Rectum")

# want subset with only biopsy samples, no follow-up measurements, no missing values in key covariates
#no.follow <- meta.sorted[!duplicated(meta.sorted$anonymized_name),]
#meta.site <- subset(no.follow, collection=="RISK")
#meta.sub <- subset(meta.site, type_sample=="biopsy" & biopsy_location=="Rectum")

# match otu table with cleaned meta data
otu.sub <- subset(otu_table_col, rownames(otu_table_col) %in% metadata$baylor_other_id)


###Prevalence filtration & decontam was not processed as phyloseq was alread filtered
# filter otu table 

pa.table <- 1*(otu.sub > 0) #Presence absence table
#otu.sub.f <- otu.sub[, -which(colSums(pa.table)<5)]  #Prevalence filtering
#dim(otu.sub)
#dim(otu.sub.f)    # 16 taxa were removed

# remove obs with < 5000 reads


# remove obs with < 5000 reads
otu.full <- otu.sub[rowSums(otu.sub)>5000,] #samples with too low reads were removed
#Subsetting samples with effective reads
meta.sub <- subset(meta.sorted, meta.sorted$baylor_other_id %in% rownames(otu.full))
dim(otu.sub)

## [1] 35 45

dim(otu.full) #14 samples were thrown out

## [1] 35 45

min(rowSums(otu.full))

## [1] 173098

meta.full <- subset(meta.sub, meta.sub$baylor_other_id %in% rownames(otu.full))
#Making data of t-tests - Not usefull for my analysis
#Anyway ran for differences between control vs treatment
p.t.full = t.test(rowSums(otu.full) ~ meta.full$control, var.equal=FALSE)$p.value
p.wilcox.full = wilcox.test(rowSums(otu.full) ~ meta.full$control)$p.value

# --------------------------
# a subsample of 50 subjects
# --------------------------

set.seed(seed)
otu.sam10 <- otu.full[sort(sample(nrow(otu.full), 10)), ]  
meta.sam10 <- subset(meta.sub,meta.sub$sample_name %in% rownames(otu.sam10))

#1st. We will check all the samples, non-rarefied

########################
# prepare data for anal
########################

if (data_full==1) {
    otu.tab = otu.full
    meta.data = meta.full
} else {
    otu.tab = otu.sam10
    meta.data = meta.sam10
}
dim(otu.tab)

## [1] 35 45

dim(meta.data)

## [1] 35 66

########################
# permanova-F
########################
#LDM-F: averaging the F statistics over all rarefactions, 


if (data_full==1) {
    formula = otu.tab ~ treatment
} else {
    formula = otu.tab ~ treatment
}

res.F <- permanovaFL(formula = formula,
                     data=meta.sub,dist.method = dist_method,
                     scale.otu.table = FALSE, n.rarefy=100, n.rej.stop=100, seed=seed,
                     perm.between.type = "none", perm.within.type = "free", cluster.id = "subject_id")

## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000

        #Permutation within requires the same sie of clusters - as some of samples have failed sequencing this option is not valid.
# subject_id was added as a fixed term

p_permanovaF <- res.F$p.permanova


# -----------------------------
# dist for permanova-D2 (1-100)
# -----------------------------

n_rarefy = 100
label = colnames(otu.full)

n_sam_new = nrow(otu.tab)
dist_sq_all = array(0, dim=c(n_sam_new, n_sam_new, n_rarefy))

set.seed(seed)

for (r in 1:n_rarefy) {
    otu_table_rarefy <- Rarefy(otu.tab)$otu.tab.rff
    otu_table_rarefy1 <- (otu_table_rarefy>0)*1
    colnames(otu_table_rarefy1) = label
    colnames(otu_table_rarefy) = label
    if (dist_type == 5 ) {
        dist <- vegdist(x=otu_table_rarefy1, method=dist_method)
    } else if (dist_type == 4) {
        dist <- GUniFrac(otu_table_rarefy1, crohn.tree, alpha=c(1))$unifrac[,,"d_UW"]
    }
    dist <- as.matrix(dist)
    
    if (r==1) {
        dist_sq_all[,,1] = dist^2
    } else {
        dist_sq_all[,,r] = dist_sq_all[,,r-1] + dist^2
    }
}

# --------------------------
# dist for permanova-D2-A
# --------------------------
#LDM-A: averaging the RSS terms over all rarefactions, 
if (dist_type == 5 ) {
    res <- jaccard.mean(otu.tab)
    dist_sq_inf_o1 = res$jac.mean.sq.o1
    dist_sq_inf_o2 = res$jac.mean.sq.o2
} else if (dist_type == 4) {
    res <- unifrac.mean(otu.tab, crohn.tree)
    dist_sq_inf_o1 = res$unifrac.mean.sq.o1
    dist_sq_inf_o2 = res$unifrac.mean.sq.o2
}
#Negative controls were having issues - 0 reads after rarefying. These OTU table with zeros caused problem in data processing, therefore the Neg.s were removed from rarefaction analysis


#########################
# permanova-D2, D2-A
#########################

if (data_full==1) {
        
        formula = dist_sq_r ~ treatment

} else {
    formula = dist_sq_r ~ treatment
}

p_permanovaD2 = NULL
p_permanovaD2A2 = NULL
p_permanovaD2A1 = NULL

for (r in c(1:10, 50, 100, 100^2, 100^3)) { # "100^2" corresponds to D2-A2, "100^3" corresponds to D2-A1
    cat("number of rarefaction:", r, "\n")

    if (r == 100^2) {
        dist_sq_r = dist_sq_inf_o2
    } else if (r == 100^3) {
        dist_sq_r = dist_sq_inf_o1
    } else {
        dist_sq_r = dist_sq_all[,,r]/r
    }
    res_permanovaD2 = permanovaFL(formula = formula, 
                                  data = meta.data, n.perm.max=5000, 
                                  square.dist=FALSE, n.rarefy=0, n.rej.stop=100, seed=seed,
                                  perm.between.type = "none", perm.within.type = "free", cluster.id = "subject_id")
    
    if (r == 100^2) {
        p_permanovaD2A2 = res_permanovaD2$p.permanova
    } else if (r == 100^3) {
        p_permanovaD2A1 = res_permanovaD2$p.permanova
    } else {
        p_permanovaD2 = c(p_permanovaD2, res_permanovaD2$p.permanova)
    }
}

## number of rarefaction: 1 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 2 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 3 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 4 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 5 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 6 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 7 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 8 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 9 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 10 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 50 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 100 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 10000 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 1e+06 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000

#tab = c(p_permanovaF, p_permanovaD2, p_permanovaD2A2, p_permanovaD2A1, i_seed)
#write.table(t(tab), file = paste("", data_full, "_dist", dist_type, ".txt", sep=""), 
#            quote=FALSE, append=TRUE, row.names=FALSE, col.names=FALSE)


p_permanovaF

##    rarefy1    rarefy2    rarefy3    rarefy4    rarefy5    rarefy6    rarefy7 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##    rarefy8    rarefy9   rarefy10   rarefy11   rarefy12   rarefy13   rarefy14 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy15   rarefy16   rarefy17   rarefy18   rarefy19   rarefy20   rarefy21 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy22   rarefy23   rarefy24   rarefy25   rarefy26   rarefy27   rarefy28 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy29   rarefy30   rarefy31   rarefy32   rarefy33   rarefy34   rarefy35 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy36   rarefy37   rarefy38   rarefy39   rarefy40   rarefy41   rarefy42 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy43   rarefy44   rarefy45   rarefy46   rarefy47   rarefy48   rarefy49 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy50   rarefy51   rarefy52   rarefy53   rarefy54   rarefy55   rarefy56 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy57   rarefy58   rarefy59   rarefy60   rarefy61   rarefy62   rarefy63 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy64   rarefy65   rarefy66   rarefy67   rarefy68   rarefy69   rarefy70 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy71   rarefy72   rarefy73   rarefy74   rarefy75   rarefy76   rarefy77 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy78   rarefy79   rarefy80   rarefy81   rarefy82   rarefy83   rarefy84 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy85   rarefy86   rarefy87   rarefy88   rarefy89   rarefy90   rarefy91 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy92   rarefy93   rarefy94   rarefy95   rarefy96   rarefy97   rarefy98 
## 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 0.00119976 
##   rarefy99  rarefy100 
## 0.00119976 0.00119976

p_permanovaD2

##  [1] 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996
##  [7] 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996

p_permanovaD2A2

## [1] 0.00019996

p_permanovaD2A1

## [1] 0.00019996

res_permanovaD2$R.squared

## [1] 0.3094375

Rarefaction stratified (Sputum)

As a sanity check, rarefaction based analysis was conducted using LDM

if(length(args)==0){
    print("No arguments supplied")
    
    dist_type <- 5      # 1:"horn"; 5: "jaccard"; 4: "unwt-unifrac"
    data_full <- 1      # 1: full sample; 0: a subsample of 50 subjects
    i_seed <- seed         # seed for generating different subsamples
    
} else{
    for(i in 1:length(args)){
        eval(parse(text=args[[i]]))
    }
}

## [1] "No arguments supplied"

if (dist_type==5) {
    dist_method = "jaccard"
} else if (dist_type==4) {
    dist_method = "unwt-unifrac"
} else if (dist_type==1) {
    dist_method = "horn"
}  


###############################################
# read in data
###############################################

# read in otu table, taxonomy info
#dat <- read_biom("rawData/otu_table.biom")
#otu_table <- as.data.frame(as.matrix(biom_data(dat)))
#taxonomy <- observation_metadata(dat)
#crohn.tree <- read.tree("rawData/insertion_tree.relabelled.tre")

otu_table_col <- as.data.frame(t(otu_table(phyloseq$phyloseq_count %>%
                                                   subset_samples(sample_type %in% c("Sputum"))))) #make samples rows, OTUs columns
otu_table_col <- otu_table_col[,colSums(otu_table_col>0)>0] # Removing data with zeros

# read in metadata
metadata <- sample_data(phyloseq$phyloseq_count %>%
                                                   subset_samples(sample_type %in% c("Sputum"))) %>% data.frame()

# create dichotomous variable for diagnosis
#metadata <- cbind(metadata,ifelse(metadata$diagnosis=="CD" | metadata$diagnosis=="IC" | metadata$diagnosis=="UC", 1, 0))
#colnames(metadata)[61] <- "disease" 

# sort metadata to match otu table
meta.sorted <- metadata[order(match(metadata$baylor_other_id, rownames(otu_table_col))),]

# remove duplicates at the end
#meta.site.dup <- subset(meta.sorted, collection=="RISK")
#meta.sub.dup <- subset(meta.site.dup, type_sample=="biopsy" & biopsy_location=="Rectum")

# want subset with only biopsy samples, no follow-up measurements, no missing values in key covariates
#no.follow <- meta.sorted[!duplicated(meta.sorted$anonymized_name),]
#meta.site <- subset(no.follow, collection=="RISK")
#meta.sub <- subset(meta.site, type_sample=="biopsy" & biopsy_location=="Rectum")

# match otu table with cleaned meta data
otu.sub <- subset(otu_table_col, rownames(otu_table_col) %in% metadata$baylor_other_id)


###Prevalence filtration & decontam was not processed as phyloseq was alread filtered
# filter otu table 

pa.table <- 1*(otu.sub > 0) #Presence absence table
#otu.sub.f <- otu.sub[, -which(colSums(pa.table)<5)]  #Prevalence filtering
#dim(otu.sub)
#dim(otu.sub.f)    # 16 taxa were removed

# remove obs with < 5000 reads


# remove obs with < 5000 reads
otu.full <- otu.sub[rowSums(otu.sub)>5000,] #samples with too low reads were removed
#Subsetting samples with effective reads
meta.sub <- subset(meta.sorted, meta.sorted$baylor_other_id %in% rownames(otu.full))
dim(otu.sub)

## [1]  30 165

dim(otu.full) #14 samples were thrown out

## [1]  30 165

min(rowSums(otu.full))

## [1] 40700

meta.full <- subset(meta.sub, meta.sub$baylor_other_id %in% rownames(otu.full))
#Making data of t-tests - Not usefull for my analysis
#Anyway ran for differences between control vs treatment
p.t.full = t.test(rowSums(otu.full) ~ meta.full$control, var.equal=FALSE)$p.value
p.wilcox.full = wilcox.test(rowSums(otu.full) ~ meta.full$control)$p.value

# --------------------------
# a subsample of 50 subjects
# --------------------------

set.seed(seed)
otu.sam10 <- otu.full[sort(sample(nrow(otu.full), 10)), ]  
meta.sam10 <- subset(meta.sub,meta.sub$sample_name %in% rownames(otu.sam10))

#1st. We will check all the samples, non-rarefied

########################
# prepare data for anal
########################

if (data_full==1) {
    otu.tab = otu.full
    meta.data = meta.full
} else {
    otu.tab = otu.sam10
    meta.data = meta.sam10
}
dim(otu.tab)

## [1]  30 165

dim(meta.data)

## [1] 30 66

########################
# permanova-F
########################
#LDM-F: averaging the F statistics over all rarefactions, 


if (data_full==1) {
    formula = otu.tab ~ treatment
} else {
    formula = otu.tab ~ sample_type + treatment
}

res.F <- permanovaFL(formula = formula,
                     data=meta.sub,dist.method = dist_method,
                     scale.otu.table = FALSE, n.rarefy=100, n.rej.stop=100, seed=seed, perm.between.type = "none",
                     perm.within.type = "free", cluster.id = "subject_id")

## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000

        #Permutation within requires the same sie of clusters - as some of samples have failed sequencing this option is not valid.
# subject_id was added as a fixed term

p_permanovaF <- res.F$p.permanova


# -----------------------------
# dist for permanova-D2 (1-100)
# -----------------------------

n_rarefy = 100
label = colnames(otu.full)

n_sam_new = nrow(otu.tab)
dist_sq_all = array(0, dim=c(n_sam_new, n_sam_new, n_rarefy))

set.seed(seed)

for (r in 1:n_rarefy) {
    otu_table_rarefy <- Rarefy(otu.tab)$otu.tab.rff
    otu_table_rarefy1 <- (otu_table_rarefy>0)*1
    colnames(otu_table_rarefy1) = label
    colnames(otu_table_rarefy) = label
    if (dist_type == 5 ) {
        dist <- vegdist(x=otu_table_rarefy1, method=dist_method)
    } else if (dist_type == 4) {
        dist <- GUniFrac(otu_table_rarefy1, crohn.tree, alpha=c(1))$unifrac[,,"d_UW"]
    }
    dist <- as.matrix(dist)
    
    if (r==1) {
        dist_sq_all[,,1] = dist^2
    } else {
        dist_sq_all[,,r] = dist_sq_all[,,r-1] + dist^2
    }
}

# --------------------------
# dist for permanova-D2-A
# --------------------------
#LDM-A: averaging the RSS terms over all rarefactions, 
if (dist_type == 5 ) {
    res <- jaccard.mean(otu.tab)
    dist_sq_inf_o1 = res$jac.mean.sq.o1
    dist_sq_inf_o2 = res$jac.mean.sq.o2
} else if (dist_type == 4) {
    res <- unifrac.mean(otu.tab, crohn.tree)
    dist_sq_inf_o1 = res$unifrac.mean.sq.o1
    dist_sq_inf_o2 = res$unifrac.mean.sq.o2
}
#Negative controls were having issues - 0 reads after rarefying. These OTU table with zeros caused problem in data processing, therefore the Neg.s were removed from rarefaction analysis


#########################
# permanova-D2, D2-A
#########################

if (data_full==1) {
        
        formula = dist_sq_r ~ treatment

} else {
    formula = dist_sq_r | sex ~ disease
}

p_permanovaD2 = NULL
p_permanovaD2A2 = NULL
p_permanovaD2A1 = NULL

for (r in c(1:10, 50, 100, 100^2, 100^3)) { # "100^2" corresponds to D2-A2, "100^3" corresponds to D2-A1
    cat("number of rarefaction:", r, "\n")

    if (r == 100^2) {
        dist_sq_r = dist_sq_inf_o2
    } else if (r == 100^3) {
        dist_sq_r = dist_sq_inf_o1
    } else {
        dist_sq_r = dist_sq_all[,,r]/r
    }
    res_permanovaD2 = permanovaFL(formula = formula, 
                                  data = meta.data, n.perm.max=5000, 
                                  square.dist=FALSE, n.rarefy=0, n.rej.stop=100, seed=seed,
                                  perm.between.type = "none", perm.within.type = "free", cluster.id = "subject_id")
    
    if (r == 100^2) {
        p_permanovaD2A2 = res_permanovaD2$p.permanova
    } else if (r == 100^3) {
        p_permanovaD2A1 = res_permanovaD2$p.permanova
    } else {
        p_permanovaD2 = c(p_permanovaD2, res_permanovaD2$p.permanova)
    }
}

## number of rarefaction: 1 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 2 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 3 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 4 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 5 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 6 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 7 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 8 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 9 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 10 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 50 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 100 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 10000 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000 
## number of rarefaction: 1e+06 
## permutations: 100 
## permutations: 200 
## permutations: 300 
## permutations: 400 
## permutations: 500 
## permutations: 600 
## permutations: 700 
## permutations: 800 
## permutations: 900 
## permutations: 1000 
## permutations: 1100 
## permutations: 1200 
## permutations: 1300 
## permutations: 1400 
## permutations: 1500 
## permutations: 1600 
## permutations: 1700 
## permutations: 1800 
## permutations: 1900 
## permutations: 2000 
## permutations: 2100 
## permutations: 2200 
## permutations: 2300 
## permutations: 2400 
## permutations: 2500 
## permutations: 2600 
## permutations: 2700 
## permutations: 2800 
## permutations: 2900 
## permutations: 3000 
## permutations: 3100 
## permutations: 3200 
## permutations: 3300 
## permutations: 3400 
## permutations: 3500 
## permutations: 3600 
## permutations: 3700 
## permutations: 3800 
## permutations: 3900 
## permutations: 4000 
## permutations: 4100 
## permutations: 4200 
## permutations: 4300 
## permutations: 4400 
## permutations: 4500 
## permutations: 4600 
## permutations: 4700 
## permutations: 4800 
## permutations: 4900 
## permutations: 5000

#tab = c(p_permanovaF, p_permanovaD2, p_permanovaD2A2, p_permanovaD2A1, i_seed)
#write.table(t(tab), file = paste("", data_full, "_dist", dist_type, ".txt", sep=""), 
#            quote=FALSE, append=TRUE, row.names=FALSE, col.names=FALSE)


p_permanovaF

##    rarefy1    rarefy2    rarefy3    rarefy4    rarefy5    rarefy6    rarefy7 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##    rarefy8    rarefy9   rarefy10   rarefy11   rarefy12   rarefy13   rarefy14 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy15   rarefy16   rarefy17   rarefy18   rarefy19   rarefy20   rarefy21 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy22   rarefy23   rarefy24   rarefy25   rarefy26   rarefy27   rarefy28 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy29   rarefy30   rarefy31   rarefy32   rarefy33   rarefy34   rarefy35 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy36   rarefy37   rarefy38   rarefy39   rarefy40   rarefy41   rarefy42 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy43   rarefy44   rarefy45   rarefy46   rarefy47   rarefy48   rarefy49 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy50   rarefy51   rarefy52   rarefy53   rarefy54   rarefy55   rarefy56 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy57   rarefy58   rarefy59   rarefy60   rarefy61   rarefy62   rarefy63 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy64   rarefy65   rarefy66   rarefy67   rarefy68   rarefy69   rarefy70 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy71   rarefy72   rarefy73   rarefy74   rarefy75   rarefy76   rarefy77 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy78   rarefy79   rarefy80   rarefy81   rarefy82   rarefy83   rarefy84 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy85   rarefy86   rarefy87   rarefy88   rarefy89   rarefy90   rarefy91 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy92   rarefy93   rarefy94   rarefy95   rarefy96   rarefy97   rarefy98 
## 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 
##   rarefy99  rarefy100 
## 0.00019996 0.00019996

p_permanovaD2

##  [1] 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996
##  [7] 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996 0.00019996

p_permanovaD2A2

## [1] 0.00019996

p_permanovaD2A1

## [1] 0.00019996

res_permanovaD2$R.squared

## [1] 0.3931708